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The self-adaptive grid code, sage 

Version 3 

ijt 4. 

Carol B. Davies and Ethiraj Venkatapathy 1 
Ames Research Center 

SUMMARY 

The previous version of the SAGE code (Version 2) is described in NASA TM-1 10350 (Davies and 
Venkatapathy, 1995). This new report on Version 3 includes all the information in the original publication 
plus all upgrades and changes to the SAGE code since that time. The most significant upgrade is the feature 
to move the outer grid boundary, either reducing or expanding the grid, based on flow features or other 
information. This ability was originally featured in a separate code called OUTBOUND developed by the 
same authors and extensively used in-house. In addition, a new reclustering algorithm has been added that 
permits very tight clustering in the boundary layer, a significant asset for grids used in viscous flow 
calculations. 

The first section of this document describes the formulation of the adaption method. The second 
section is presented in the form of a user guide that explains the input and execution of the code, while the 
third section provides many examples. A special case is the ability to handle the PLOT3D BLANK option 
that is used in complex multigrid structures to blank out regions of overlapping grids. Due to its complexity 
and size, this option is not available in SAGEv3 but is described in this document in Section 5 and the 
special SAGE version can be obtained from the authors. 

Successful application of the SAGE code in both two and three dimensions for the solution of 
various flow problems has proven the code to be robust, portable, and simple to use. Although the basic 
formulation follows the method of Nakahashi and Deiwert (1985), many modifications have been made to 
facilitate the use of the self-adaptive grid method for complex grid structures. Modifications to the method 
and the simple but extensive input options make this a flexible and user-friendly code. The SAGE code can 
accommodate tw ? o-dimensional and three-dimensional, finite-difference and finite-volume, single grid and 
zonal-matching multiple grid flow problems. 


1. METHOD 
1.1 Introduction 

Solution-adaptive grid methods have become useful tools for efficient and accurate flow p predictions. 
In supersonic and hypersonic flows, strong gradient regions such as shocks, contact discontinuities, and 
shear layers require careful distribution of grid points to minimize grid error and thus produce accurate flow- 
field predictions. Frequently, the generation of the first grid topology does not adequately capture these 
flow structures. It has been shown that an effective way of obtaining accurate solutions is by intelligently 
redistributing (i.e., adapting) the original grid points based on the initial flow-field solution and then 
computing a new solution using the adapted grid. The cost efficiency of grid adaptions in terms of CPU 
time depends on the basic formulation of the adaptive-grid solver. A short historical review and a list of 
references on the adaptive grid procedure used in the SAGE code are given in Davies and Venkatapathy 
(1991). 

The self-adaptive grid procedure outlined by Nakahashi and Deiwert (1985) has evolved into a 
flexible tool for adapting and restructuring both two-dimensional (2-D) and three-dimensional (3-D) grids. 
The adaptive-grid method is based on grid-point redistribution through local error minimization. The 
procedure is analogous to applying tension and torsion spring forces proportional to the local flow gradient 
at every grid point and finding the equilibrium position of the resulting system of grid points. The 
multidimensional problem of grid adaption is split into a series of one-dimensional (I-D) problems along 
the computational coordinate lines. The reduced 1-D problem then requires a tridiagonal solver to find the 
location of grid points along a given coordinate line. Multidirectional adaption is achieved by the sequential 
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application of the method in each coordinate direction. This approach produces an extremely CPU- 
efficient algorithm. 

The tension forces direct the redistribution of points to the strong gradient regions. The torsion 
forces relate information between the families of lines adjacent to one another, to maintain smoothness 
and a measure of orthogonality of grid lines. These smoothness and orthogonality constraints are direction- 
dependent, since they relate only the coordinate lines that are being adapted to the neighboring lines that 
have already been adapted. This implies that the solutions are non-unique and depend on the order and 
direction of adaption. 

The Self-Adaptive Grid codE (SAGE) code has been built with many flexible elements that make it 
user-friendly. The second section of this report is a user guide, which is independent of the first section and 
can be used as a separate document. However, the nomenclature and details of the grid adaption procedure 
within the code consistently follow the analysis, allowing the user to understand the code and implement 
any individual changes, if desired. The user guide describes the input parameters and their effects as well as 
the adaption procedure, execution commands, and subroutines, and provides many examples. 

1.2 Formulation of Adaptive-Grid Scheme 

As stated in the introduction, adaption takes place as a series of one-dimensional adaptions. Figure 
1 illustrates this concept: three constant k planes of an initial grid are shown in Fig. 1(a) and a flow-field 
solution has been obtained using this unadapted grid. The points in this grid are then adapted to the flow 
solution, starting on the first line Q—l) on the lower plane k=l . In Fig. 1(b), the first plane has already been 
adapted and the second plane is the current adaption plane. The current adaption line (j) is shown, with 
previous lines already adapted and subsequent lines awaiting adaption. The third plane is still in its original 
form. Adaption is performed in this line-by-line, plane-by-plane manner until all requested planes are 
complete. It is then possible to perform an adaption in a second direction, adapting on top of the already 

adapted grid. The number and 
order of adaptions are arbitrary 
and depend on the type of flow 
problem and the purpose of the 
adaption. 

Figure 2 shows a 
segment of the current 
adaption line in more detail. 
The lower plane has already 
been adapted and the upper 
plane is currently being 
adapted. Four forces control 
the redistribution of points 
along each line: the two 

tension springs that act on 
each side of a node, and the 
two torsion springs that 
control the smoothness of the 
grid. The tension forces co 
cluster the redistributed points into the high-gradient regions. The torsion forces (r and vj/) restrain the 
tension forces to maintain continuity between sequentially adapted lines. As shown in Fig. 2, the x force 
acts from the previously adapted line within the current plane and the y force acts from the previously 
adapted line in the preceding plane. 

A 3-D grid is described in terms of its computational coordinate directions (i,j,k) and its physical 
coordinates (x,y,z). An adaption can take place in any or all of the computational coordinate directions and 
the SAGE code permits any combination. However, for clarity, this analysis assumes that all adaptions are 
performed in the i direction and that stepping (within a plane surface) occurs in they direction. Thus k is 
the marching direction for plane stepping. The current adaption line is thus j in the plane k where lines j-1, 
j-2, etc., have already been adapted. In addition, this report uses the term ‘plane’ to mean any 3-D 
coordinate surface. 



Figure 1. 3-D adaption, (a) Initial grid; (b) directional adaption. 
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The First step in the formulation of the adaption algorithm is to consider the adaption of a single 
line, with no torsional constraints (i.e., no smoothness control). Figure 3 shows a 2-D surface, where the 
arc length at A (i.e., with s i JJt = .0) along the current adaption line,/, is defined as 


5. = S. . + 
l l ~ 1 


(X. 


,r 


(i) 


marching direction within plane 


CO 


T3 



(extension spring between each node 
t=torsion spring within current plane 
(p=torsion spring between planes 


Figure 2. Line-by-line adaption showing tension 
and torsion forces. 


Eq. (2) produces 

X A - v / = s ' m , l v = ^X 


A tension force, is defined to act 

between each / and /+/ node such that 

co, A'', = K (2) 

where oo, the tension constant, is a weighting 
function based on the flow gradient; As, = s l+! -s ( ; 
and K is the resultant force. To redistribute the 
points along a line with the minimum solution error, 
the adaptive procedure defines the weighting 
function as proportional to the derivative of the 
flow variable (Nakahashi and Deiwert 1985). In this 
formulation, CO is defined as a function of the 
normalized flow gradient, /, such that 

co-1 + Af t * (3) 

where A and B are constants directly related to the 
grid spacing and are chosen to maintain the grid 
intervals to within the limits (A s KiJN and As max ) set 
by the user on input. A and B are defined in 
Appendix I of this section. 

Equation (2) is written for each node on the 
line, giving a 1-D formulation that can be solved 
directly for As,. Taking the sum of both sides of 


/=] 




K = s„. n ► — 


giving 


Substituting back in Eq. (2), we obtain 
/ 

i 



As. = s„ 


co 


,1 


to «/ 


(4) 


(5) 



In the SAGE code, th is l-D solution 
technique is used along the initial adaption 
line of the initial adaption plane, where no 
directional information is available. 

Continuing this approach for successive line- 
by-Iine adaptions will create a mesh that is 
insufficiently smooth for input into 
computational flow-field codes. Therefore, 
to ensure a more reasonable grid, the 
redistribution of points (driven by the 
tension springs) is constrained by torsion 

terms defined on two adjacent adapted lines: one on the current plane and one on the preceding plane. 
Within the current plane, the torsion parameter x represents the magnitude of the torsion force that 
maintains smoothness and orthogonality between the node (if k) and the nodes (ij~l,k) and (ij-2,k). This 
is the only torsion parameter available on the initial surface. For subsequent surfaces, the torsion parameter 


Figure 3. Current adaption line j on a 2-D surface. 
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\j/ constrains the movement between the same node (ij,k) and the nodes (i,j,k—I) and (iJ,k—2) on adjacent 
computational planes. The torsion terms are evaluated as x t (s' — s t ) and ~ s >) ( see Fig. 3: s is at 

A ' ), where the torsion parameters land \|f define the magnitude of the torsion forces and s'and s define 

their inclination (i.e., orthogonality and smoothness). 

To introduce the torsion forces to the system of equations, we first rewrite Eq. (2) to represent the 
force balance, 

CO ,(s M -s t )~ CO, , (s,- s,_ , ) = 0 (6) 

and then add the torsion terms to obtain 


CO ,(s M -s l )-<n l _ l (s l -s l _ x )+x l (s' i -s l )+y l (s , l -s l ) = 0 (7) 

which is rearranged to produce the coefficients of the tridiagonal matrix used to solve for , 

+ +x, + \|/Ji,+co,s, +l =-T,s'-\y l s* (8) 

This equation is written for each interval along the adaption line, producing a system of (n,-l) equations. 
Since .v, and s n are known, there are n, - 2 rows in the matrix. This equation is solved iteratively (updating 

to, at new .v,) until V. 1 1 —s^ n 1 ! <10 ^ This system of equations is central to the adaption 

technique used in the SAGE code and most of the description that follows pertains to deriving the 
coefficients of this equation. 


1.3 Auxiliary Equations 

1.3.1 Tension parameter, co f . As described in the formulation, to, acts as a tension force in the 
interval (s l+i —s,J and can be imagined to be a spring (aligned with the grid line) connecting the two grid 
points, as shown in Fig. 3. This tension parameter is defined in Eq. (3) as 

co, = 1 + Af* 

where f, is a function of the gradient of the flow variable, q : 

J f- fnm_ and j = jjL = ~ (9) 

J f - f ds As, 

J max J mm 1 

The standard format of the user’s input flow-field file contains the conservative flow variables Q, as 
defined in the plotting software package PLOT3D (Walatka et al., 1990). This format (for a 3-D dataset) 
assumes the flow variables are p, pu, pv, pw, and e. Since the adaption is based on a scalar function q , this 
function can be evaluated as a user- specified combination of flow variables Q. Pressure, Mach number and 
temperature ratio are also provided as adaption variables and are internally computed, assuming the perfect 
gas relationship. However, Q need not be restricted to conservative flow variables, but can be any 
combination of user-specified flow variables that represent the flow field. 

To remove the unwanted oscillations from the discrete evaluation, f t is smoothed by adding a 
second derivative term, i.e., f, =.75/ + .125(/) + , By default, two smoothing passes are made, but 

this can be overridden by changing the filter input control parameter, NFrLT. The tension parameter, co,, 
is smoothed in the same way. 

1.3.1. 1 A and B. As seen in Eq. (3), co is also a function of A and B. These variables are the 
important elements of the self-adaptive nature of the scheme because they control the size limits of the 
grid spacing. A and B are computed from the user-supplied maximum and minimum allowable grid spacings, 
A s UAY and Av M/(V . These are relative values of mesh size and allow the user to specify the proportion of 
largest to smallest As in the final adapted grid. 

The value of A is constant throughout the grid adaption and is given by 

A = te*WL-\ ( 10 ) 

^ MIN 
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The value of B is computed (by an iterative process) for each j line to provide 

input As = computed As „„„ 

That is, the computed minimum grid spacing is equal to the user-requested value. Appendix I gives a detailed 
description of the derivation of these two parameters. 

1.3. 1.2 Torsional effect on the evaluation of CO. As just noted, the function of 
variables A and B is to control the size limits of the adapted-grid mesh spacing and, as described in 
Appendix I, B is computed from the 1-D Eq. (2). Thus the addition of the torsion terms is somewhat 
inconsistent with the value of B , and the resulting grid spacings may give values slightly outside the 
requested limits. To provide for additional control, a weighting factor is applied to co such that 

co , ■,= (1 + 4// H°/ * where co, is the weighted tightness factor. When equation (8) is solved, each As, is 
tested and if it does not lie within the requested spacing limit, co, will be computed; otherwise co, =1.0. If 
Av, < As WN , we need to relax (decrease) the value of co, to produce a larger As, in the next iteration. If As, 
is too large, co, should be increased. We therefore use the modifier 


1 , As, 1 As, 

co, ! -+l)or co, =-(— J — + \) 


2 As 


MIN 


2 As 


MAX 


depending on whether As,<As Km or As, > As MAX . Equation (8) is therefore solved using the modified 
values of co. It should be noted that co, =1.0 in the regions of side-spacing control (see section 1.4.3) 
where As, is permitted to be outside the specified limits. 


1.3.2 Torsion terms. As shown in Eq. (7), the torsion terms, x/s'-s,) and \\f f (s*-s,)> are added 
to the 1-D equation to preserve the required smoothness and orthogonality of the grid. They are 
constructed as being analogous to the behavior of torsion springs. Consider first the torsion term within the 
plane, / (x). A torsion force, T, relates the node at (i,j-l,k) to the node (ij,k), is proportional to the 
turning angle, 0, shown in Fig. 3, and is defined as 

T, = ( 11 ) 

where K is a torsion-related constant. The next step is to relate this torsional force to the previously 
defined parameters on the adaption line. Since 0 is generally small, we can approximate 0, by 
(s' -s f )j\DA'\ (see Fig. 3), where (s' — Sj) is a function of the inclination angle computed from the 

proportioned orthogonality and straightness vectors, as described in section 1.3. 2. 3. In addition, we assume 
that K is proportional to the maximum co, (= I + A) and the local aspect ratio of the grid cell. Hence we can 
write 


— max ( S t+\,j-\ \) 

K ” i\da'\ 


( 12 ) 


where X is the proportionality coefficient and is an input quantity. The variable x used in Eq. (8) can thus 
be defined as 


T: =’ 




J - » 




(13) 


2 \DA'\ 


The evaluation of \j f is similar, with the node (i f j~I,k) replaced by (ij,k-l) and a proportionality 
coefficient X defined that is analogous to X such that: 


V, - 


^ ® max /+ j j jf . 


1 *i-\j,k-\ 


) 


(14) 


2 \BA'\ 
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where BA' is the distance from (ij,k-l) to the location of s* j k . The evaluations of \DA'\ and BA' are 


given in Appendix II. 



1.3.2. 1 Evaluation of s' and s*. These two variables, seen on the right-hand side of 
Eq. (8), are calculated from the torsion vectors and provide the direction of the smoothness and 
orthogonality constraints. The location of s' lies at the intersection of the projection of the within-plane 
torsion vector 7, with the arc-length vector s, along the current adaption line. The location of s * is 
similarly determined from the between-plane torsion vector t* . Figure 4 shows the two unit torsion 
vectors: 7, acting from point D and intersecting the line AC at A', and (*, acting from point B and 
intersecting the same line at A *, such that 

s' = s, + AA' and s* = s, + AA* (15) 


In general, neither 7 nor 7* will directly intersect the line AC, and the projection of the torsion vectors 
onto s is required. An example of the projection vectors is shown in Fig. 5 where the projection of t onto 
the k plane intersects s at A' . A full description of the calculation of the intersection is given in Appendix 
II. 

1.3.2. 2 Torsion vectors t and t*. Figure 6 shows in detail the within-plane torsion 
vector 7, and the associated base vectors (n = f(u,b), and e). The vectors in these figures are displayed on 
a 2-D surface, and it should be remembered that the shown 7 is actually the projection of 7 seen in Fig. 5. 
The unit vector 7, associated with the nodal point (ij,k) [but acting from the point (i,j-I,k)] is defined as 

Kj.k= 06) 

where 


C is the input parameter that defines the percentage of straightness to orthogonality, 
e is an average straightness vector from ( i,j-2,k ) to ( i,j-J,k ), and 

represents the orthogonality vector between the j-1 line and the node (ij,k) and is a function of 


b and it vectors described below. 

The between-plane torsion vector 7* which acts from the node (i,j,k-l), is similarly defined as 




, =A/c;*;+(i-c>;/ 


(17) 




where C* is the input proportion parameter, and 
e* and n* are the equivalent straightness and 
normal vectors. 

It should be noted that Figs. 4 and 6 show 
vectors intersecting the j line in the interval 
(i,i+l)\ however, this is not necessarily the case. 
In general, we assume that the projection of t i 
intersects the j line in the interval (1,1+1). 
Appendix II describes a general method used to 
compute the intersection of a vector with the 
vector s. 

1.3. 2.3 Evaluation of vectors s, e, and n. 
The unit vector s i (direction AC in Fig. 7) is the 
arc-length vector from s i to s j+[ along a j line. It 
is defined as 

hj= s .J + s y,J + S z) 

where 

^ (^i+xjjc Sy (yi+\j,k y,j,k) 

s, s, 



Figure 6. Torsion vector as a combination of 
normal and straightness vectors. 


and =n 

hi 

The unit vector e i (direction ED in Fig. 7) is defined as the sum of three straightness vectors, 
d x _ j, d t , and ^ where 

Kj =d x,‘ +d y,j + d J 


such that 



Figure 7. Straightness vectors, e and e * . 


=^( X i.j-U- X '.j-2.k)’ d y, =^j\( y i.j-l.k-yi.J-2.k) 


and d *' = 17] r Z ‘J- ] k ~ z ' J- 2 - k ) 

Therefore, e t r +d t +d l+l ). If the line 

V3 

j-2 does not exist (such as at the second line 
from a physical boundary), it is assumed that 
e — h . 

The straightness vector between planes, 


e% is shown as PB in Fig. 7. It is computed like 
e above, with points (ij-l,k) and ( ij-2,k ) 
replaced with (ij,k-l) and (if k-2) respectively. 

The vector h is a combination of two 


vectors: 



+ (\-t n )u 


( 18 ) 


where u i represents orthogonality to the ik plane through the j-1 line, and b f the orthogonality to the Ik 

plane through the j line, as shown in Fig. 8. The parameter t n that proportions b and it is defaulted to .5, 
changing only at the upper and lower marching boundaries, as discussed in section 1.4.4. The calculation of 
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a normal vector to a point in a 3-D 
surface needs defining: for example, in 
Fig. 8, the vector u acts at D , 
perpendicular to the ik surface at j—1 (the 
surface that contains G, F f M, and Q). In 
this case, there are four possible planes to 
choose from {MDF, MDQ, GDF \ and 
GDQ ) and thus the definition of the u 
and b vectors is the average of the 
normals to each of the four planes, when 
the planes exist. Since two vectors define 
the plane in which they lie, the normal 
to that plane is obtained by taking their 
cross product. The vector b would also 
act at D but is perpendicular to the ik 
plane at j. Appendix III describes the 
general calculation of the normal to a 
given plane at a given point. 

Tn addition to the within-plane 
orthogonal constraints, orthogonality is 
also controlled between planes by a 

function of the vector n* , defined as 




(19) 


The normal u* acts at point B, perpendicular to the ij plane at k-I, and b* acts at A , but perpendicular to 
the plane at k . 


1,4 Treatment of Boundaries 

The ability to modify the adaption techniques in boundary regions substantially improves the 
flexibility of the adaptive scheme. The adaption domain is defined by the user and may be equal to, or a 
subset of, the physical grid. A boundary occurs at the limits of the adaption domain defined by the user. 
Within a plane, there are two types of boundaries: marching boundaries (all points along the initial and final 


adaption lines) and edge boundaries (at i — i st and 
i = i efui ). There are also initial and final surface 

boundaries to take into consideration. Figure 9 shows 
a 2-D example of an adaption domain as a subset of 
the physical grid, and illustrates the two types of 
boundary lines w hen marching in the j direction. Note 
that if marching had been in the / direction, the 
boundary definitions would have been reversed. 

By internally amending the previously defined 
variables of C r t n , and X in the boundary regions, it 

is possible to maintain physical characteristics at 
boundaries, allows for multiple passes, provide 
continuity across grid boundaries, and simplify 
multiple grid adaptions. Specialized boundary 
conditions, such as wall shape preservation and 
periodic boundaries, are discussed in separate sections. 

1.4.1 Treatment of initial marching 
boundary line. If the initial adaption line is within 
the physical domain, a smooth transition will be 



physical domain. 
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required across the starting internal boundary. For example, this situation occurs when adapting in zones; 
each zone has different flow features and the user may wish to march up to a certain line using one set of 
parameters, then continue marching using a new set. The common boundary between the two zones must 
remain unchanged when starting the second adaption pass. This feature is controlled by the input parameter 
m g . When m g >0, a smooth transition from external grid lines (e.g., those in the already adapted zone) to 


internal lines is created by maintaining the same grid points along the initial line and then incrementally 
introducing the input adaption parameters. For example, when stepping to the second adaption line, most 
straightness is maintained by setting C, = I (see Eq. (16)). At each subsequent line, C f is gradually decreased 
until it equals the user-supplied input value. The number of lines stepped before the full effect of the new 
parameters is felt depends on m g \ C t will linearly decrease until m lines have been adapted. An example of 

this can be seen in Fig. 9, in which j st =4 and m g =5. The grid points along j = 4 will remain unchanged 


while those along j = 8 will be fully adapted to the input control parameters. The code controls the 
adaption parameters for lines in between. Consequently, we define a variable n m as 


n 


m 


Hl n 




■Js, 


( 20 ) 


and replace the value of C, used in Eq. (16) by 


C, 




At the same time, the value of A. is increased to A(1 +5 n m ) so that this amendment to C t is more effective. 
After m g lines, n m = 0, thus returning C { and A to their original values. Note that the computation of e 
along the / vf + 1 line is a function of the j st — 1 line (if it exists), and this will also help in the merging 
process. 

A similar parameter, m*, is used when the first plane must remain unchanged (e.g., at a multigrid 
boundary) and subsequent planes are gradually adapted. 

1.4.2 Treatment of final marching boundary line. Another discontinuity will occur between 
the final adaption line and any subsequent lines external to the adaption domain. Since the adaption process 
is a marching scheme, it is not possible to use the same merging concept described earlier. On request 
(through the MARCH input parameter), the grid points on the remaining external lines will be redistributed 
with the same proportions as the points on the final adapted line. In addition, nonadapted planes can be 
proportioned with respect to lines on previously adapted planes. This is not an adaption to the flow field, 
but provides a more acceptable interface between the computational and physical domains. 

1.4.3 Treatment of side-edge boundary. A smooth transition will also be needed at the side- 
edge boundaries if the adaption domain is internal to the given grid, i.e., if i st > 1 or i t , mJ < i max . Figure 
10(a) shows an example in which the 
fixed external grid spacing is denser than 
the first redistributed points, giving a 
discontinuous effect across the boundary. 

If the user requests continuity of mesh 
spacing across the side-edge boundaries 
(by setting the input parameter NEDGE, 

(n g ) ^ 0), a modification is made to the 

tension parameter, co. Consider the start- 
edge boundary' at / - / v/ along line j. We 
wish to enforce some value on co, that 
will give a value of As, close to the value 
of To do this, we find the average 

cdAv along the converged j-1 line and 
replace co, by coA.y/As, , } This value remains fixed during the iteration, but is merged into the updated 
values of co close to the boundary. In general, this implies that co 2 , co 3 , and co 4 will be amended; however, 
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Figure 10. Control of side-edge boundaries. 

(a) With no edge control; (b) with edge control 
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an additional option is available (using input parameters MG1 and MG2) to spread the effect further into 
the adaption domain. Figure 10(b) shows the effect of this process, giving a more appropriate spacing in 
the vicinity of the i s , boundary. The end-edge boundary is handled in a similar manner. 

The side-edge spacing often needs to be improved even when the adaption domain coincides with 
the physical domain. If an adaption pass generates inappropriate side-edge spacing, the code can be rerun 
with n K set to nonzero in an attempt to improve the spacing by using the above technique. In this case we 

do not have an external As, and A s 2 j is used instead of As,- a _|j. This will usually prevent the spring 

constants from pulling the lines too far off the boundary. 

The variable can be set to initiate computation for either or both side edges. 

1.4.4 Treatment of orthogonality at boundaries. The code provides the choice of 
constructing grid lines that are as orthogonal as possible to a marching boundary, either to the final line in a 
plane or to the entire final plane. To accomplish this, the normal vector h (and/or «*) is emphasized over 
the straightness vector by decreasing C, when the adaption line is close to a marching boundary line. For 

even greater control, the coefficient t n is modified to emphasize either u or b, depending on whether the 
initial-line or end-line boundary is being considered. To ensure that the modifications to these torsion 
coefficients sufficiently affect the computation, the value of X is simultaneously increased to accentuate 
the torsion term. Since not all marching boundaries are physical boundaries, an input option is available 
that will override this emphasis on orthogonality for either boundary. 

1.5 Preservation of Initial Wall Boundary Shape 

In applications for which the shape of the wall boundary is defined by large geometric gradients, 
such as for turbine blades, sharp corners, and the leading edges of airfoils and wings, sufficient points need to 
be placed in the appropriate regions of the initial grid to accurately define the geometry. If flow-field 
gradients are weak in these regions, the standard grid-redistribution algorithm will cause points to be 
dispersed, leaving insufficient points to properly describe the surface. To maintain necessary clustering in 
these regions, on user request a new variable is introduced that is a function of the geometry gradients that 
define .the boundary shape. The weighting parameter now becomes a function of both flow-field gradients 
and geometry gradients. The solution procedure will therefore redistribute the points into regions of high 
geometry gradients as well as into regions of high flow-field gradients. Both the start and end wall 
boundaries are treated. A merging technique is used to integrate the boundary and internal redistribution, so 
that the internal flow is controlled only by the flow-field gradients. 

The original weighting function was defined in Eq. (3) as (0 = 1 + Af B where / is a function of 
dq/ds. When the geometry option is invoked, f becomes a function of both dq/ds and dg/ds, where 
f (dg/ds) is equal to the radius of curvature R at the wall boundaries. The radius of curvature is defined as 


where the derivatives are computed from the spline coefficients. 

The full effect of the geometry function should be felt at the wall boundaries; it should not be a 
factor in the internal grid redistribution. Flence the final weighting function takes the form 

co = 1 + Af B 

where / - C l/ f(q) + C^f(g) and the constant C (/ normally equals I. The constant C g equals I when the 
upper and lower wall boundaries are being adapted, but it must be decreased away from the walls until C g = 0 
internally. The value of the aspect ratio was chosen to drive the rate of decrease of C g \ i.e., C g = 1-4A V , 
where A s is the aspect ratio at the maximum radius of curvature. Regardless of the value of A s , C g remains 

positive for a minimum of four steps from the boundaries and decreases smoothly. When the upper 
boundary is approached, C must be turned on when necessary and increased with each step, and f(g) 

becomes a function of the upper wall boundary only. 


dx fVy dy 3 2 .v^ 
ds ds 2 3.v ds 2 


/? = 
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It is also possible to adapt only to the geometry gradient. If this is requested, C q = 0 and C g = 1 
throughout the flow, and f(g) is computed for each adaption line, based on the local geometry of that line. 
This option has proven to be a useful tool for improving the starting grid before any solution is obtained. 
Points will be smoothly clustered with respect to the geometry gradients. 

1.6 Finite-Volume Grids 

In applications using finite-volume techniques, each flow-field variable ( q ) is evaluated at the cell 
center instead of the nodal point of the grid. These cell centers are usually positioned at the average of the 
surrounding grid points (four for 2-D and eight for 3-D). This creates a solution file that has one less point 
in each direction than the grid file. The finite volume module in SAGE interpolates the q values onto the 
grid points before performing the adaption. On conclusion the q values are re- interpolated at the new cell 
centers. 

A feature of finite volume grids is the definition of the outermost boundary cells in all 
computational directions as ‘ghost cells’ containing specified flow values. The adaption procedure destroys 
these values and SAGE employs a very simple method to replace them: on output, each ghost cell contains 
the same value of q as the adjacent internal cell. It is therefore highly recommended that the user re-apply 
the boufidary conditions before processing the new grid. This is true for both 2-D and 3-D grids and implies 
that if a user wishes to adapt one plane for testing purposes, an internal plane would be more appropriate. 

1.7 Multiple Grids 

When investigating the flow around complex structures, computational grids are frequently 
organized in multiple-grid format. This enables each individual grid to be of manageable size while 
maintaining a single dataset and also allows for overlapping grids. The original SAGE code could only adapt 
single grids and it was therefore necessary to separate the grids before adaption. The latest version of SAGE 
can read and write datasets in multigrid format and provides an input control parameter to specify which 
grid to adapt. In addition, a feature is available to transfer data between matching zonal boundaries (with 1 : 1 
mapping) in separate grids. Some or all boundaries of each individual grid will match in some part to 
boundaries in one or more of the other grids. It is important that common boundaries, where data in 
separate grids represent the same location in computational space, retain the same grid and flow 
distributions. This can be seen in the simple two-part multiple grid shown in Fig. 1 1(a), where the shaded 
plane is common to both grids. After adaption, both planes must still contain matching data. 



Figure 11. Simple multi-zone grid, (a) Common plane; (b) restricted adaption direction with 
original SAGE; (c) preferred adaption direction , using plane transfer procedure. 

The adaption technique used in SAGE is a marching scheme, and therefore the order and direction 
of adaption have a marked effect on the final grid-point redistribution. One way to adapt the two grids and 
still retain the common boundary data is shown in Fig. 1 1(b): the common plane is adapted first, with plane 
marching occurring in opposite directions. The export/import feature enables the order of adaption shown 
in Fig. 1 1(c): after grid 1 has been adapted, adapted data from the final plane can be transferred to the first 
plane in grid 2. Then, using the merging feature (MGPLS) to prevent the first plane of grid 2 from being 
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corrupted, adaption can continue in the same plane marching direction. All three steps (adaption of grid 1, 
data transfer, adaption of grid 2) can be accomplished in the same execution pass. Figure 1 1(c) shows only 
the simplest multigrid structure. In reality, multiple grids may contain many grids, with surfaces, or subsets 
of surfaces, matching to more than one grid. Input options are available to handle this and examples can be 
found in Section 3. 


1.8 Outer boundary movement 


1.8.1 Moving the outer boundary. Once a solution has been obtained on an initial grid, it may 
be apparent that the location of the outer boundary of the grid is inappropriate. For example, a shock wave 
may have developed that passes out of the boundary; or alternatively, no features can be seen in part of the 
outer region and much of the grid is wasted. The ability to move the outer boundary is a very useful feature 
and is invoked by adding the input option MVBOUND to the input parameter list. The new boundary 
location moves along each existing grid line. If the boundary moves out, the new location lies on the 
extension of the vector joining the last two existing points. This feature has been extensively used for 
applications that compute flow solutions on a vehicle configuration flying a trajectory over a wide range of 
angles of attack. 

The user has three options for determining the new boundary. The first is the most simple: moving 
the boundary 7 point a certain percentage from its current location, either in or out. Along each adaption 
line, the value of s max is increased or decreased by a user-supplied percentage (controlled by 


MVBOUND =±/?, where p is the percentage), and then transformed back to (x,y,z) to form a new boundary 
point. 

The second option permits the user to parallel the outer shock (MVBOUND =999 . ) defined by the 
requested variable, INDQ. The outer boundary will be moved in or out for each adaption line, depending on 
the location of the shock. A search is made, starting at the outer boundary point of each grid line, and 
moving towards the inner boundary. As soon as a significant change in gradient is encountered [given as 


f(q t )- fall _,)> ■ 0001 where fiq. ) = 


], it is assumed that the shock location, s iKhlK . k , has been 


/urn- 


found. The relative distance of the new boundary point from the shock location is a function of d m , also 
provided by the user, such that s = s ishm . , k + s * d sn /1 00.0 . If a shock is moving out from the current 

boundary, this method Will only expand the grid by a given amount and may not project the shock location 
sufficiently outside the boundary . If it is obvious that the shock will still not be contained within the new 
boundary, moving the boundary out with option 1 and then applying this second option, will be more 
appropriate. Note that when the new location is dependent on a gradient location, the outer boundary may 
not be smooth and the smoothing parameters NSM and NSM1 may be required. The discontinuities occur 
since the outer regions of a grid (where a shock may be found) are often the most sparse. 

The final option is to move the boundary location based on a constant specified value of a flow 
variable, for example, a line parallel to a certain density contour. Note that only variables given in the 
input Q file can be used to move a boundary, not any internally computed values. 


1.8.2 Point redistribution after boundary movement. Once the new boundary edge has been 
located and the new (x,y,z) at s max computed, the internal grid points must be redistributed. Adaption is not 
an option: the code will automatically set NOUP=.true. A smoother point distribution is first required, then 
adaption can take place as a subsequent step. 

There are two options for reclustering points without adaption to the flow. With no other 
information, the code will automatically redistribute the points proportional to the original grid-point 
distribution, i.e., 

S i nr » ~~ ^ S nmx M 


The second reclustering method is based on the stretching algorithm of Vinokur (1983) and is invoked by 
inputting RECLUST>0 and providing information for the first and last mesh spacing along the reclustering 
line. 
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1.9 Vinokur Reclustering 


This method is appropriate for redistributing grid points when very small mesh spacings are required 
at the wall as, for example, in viscous grids. The user can specify the first spacing off the wall and this 
algorithm will create a smooth mesh-size expansion through the boundary layer. The method can be used to 
redistribute points after a boundary movement, or as a stand-alone alternative to the basic grid adaption 
procedure. 

The reclustering method is based on the 1-D stretching algorithm of Vinokur (1 983) and is invoked 
by inputting RECLUST>0 and providing information for the first and last mesh spacing along the 
reclustering line. The redistribution can apply to the entire domain (RECLUST=7«/ax) or to a portion of 
the domain (e.g., RECLUST=« will redistribute the first n points from the wall). 

The Vinokur algorithm redistributes points as = f(d K ,d e ) where d w and d e are user input 

variables that help to define the first and last As spacings needed by the algorithm. The last As spacing, 
As rmax , is calculated as 

^v m ar= d eS m ,,AnptS-\) 

where the default value of d e is 5. 

There are three options to assign a value for the first wall spacing, A s wall : 

1) default to the original As,; 

2) input d w as the required value of A s waH \ or 

3) input d w =999 and then input the first (/=/) and last ( j=jmax ) wall spacings. 
Option 3 is provided to overcome a 1-D limitation of the algorithm: that there is no smoothness 

from line to line. This is emphasized in regions where the value of A s mM is too small with respect to s m(lx . 

If the grid is expanding in size downstream, the A s miU that is appropriate near the vehicle nose may be too 
small at the outflow plane. Option 3 is provided to vary the first spacing off the wall while marching along 
the body surface. If d K =999, then d w and d must be input. These are the wall spacings at the first line 

0=1), and last line Q-jmax ) and enables d w to be computed internally for each j. If p i is the distance along 
the wall (at /=/) from the current j line to the j=l line, then 

As, . = « + bjpj where a = d Wy and b = ( d Wi - d Wf )j ^jp jmax ) 

As mentioned above, it is possible to recluster the grid points without adapting and without moving 
the outer boundary. An input of RECLUST=ft, n>3 and MVBOUND=0 will invoke the re-clustering 
algorithm between points s { and s n . Here, n could indicate the edge of the boundary layer, or the entire 
grid. If n<imax , the ‘outer’ edge spacing sent to the reclustering algorithm is - (s n+x -s n _ x )/ 2 . 

1.10 Optional Versions of SAGE 

1.10.1 Blanked grids. Certain types of complex multiple-grid structures utilize the facility in 
PLOT3D that blanks out regions of overlapping grids. However, the use of multiple, overset grids with 
blanking creates a new environment for the adaption algorithm. The adaption procedure used in the SAGE 
code is a 1-D approach; i.e., a single line is adapted before stepping to the adjacent line. The integrity of 
the adapted grid is maintained by defining restraining forces between the current adaption line and 
previously adapted lines, either within the current plane or on the previously adapted adjacent plane. The 
method assumes that there are an equal number of grid points along each line; i.e., there is a one-to-one 
correspondence between a point and its associated restraining forces. 

Overset grids are created in the normal structured manner, with equal points along each line. 
However, a flag associated with each grid point indicates whether the point should be included as an active 
part of the grid. This method creates adjacent grid lines with different numbers of active grid points, and 
even more problematical, regions of grid holes (i.e., one or more gaps along the line). The challenge for the 
SAGE algorithm and code is to develop a method that can accommodate these complications. Section 5 of 
this document contains a description of the new method that was developed and it can be seen that it adds 
complexity to the code. As a result, it was decided not to include this feature in SAGEv3 at this time. 
However, SAGEv2B (Version 2 with Blanking) is available from the authors. The description in Section 5 
applies to this code but will remain appropriate if the feature is added to a later version of SAGE. 
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1.10.2 Cyclical and periodic boundaries. Certain classes of grids used by flow-field solvers are 
not suitable for adaption using the basic grid formulation described thus far. The self-adaptive grid procedure 
is a marching scheme; i.e., the solution along each line is influenced only by previously adapted lines. The 
adaption of the first grid line is based on flow gradients ORly, but as marching proceeds, the redistribution of 
points is dampened by the torsion effect. As a result, the final adaption line will be somewhat less adapted 
to the flow-field gradient than the initial line. This implies that the scheme cannot be directly applied to 
cyclical and periodic grid structures, since common and/or matching boundaries will show as discontinuities. 
For 2-D problems, these difficulties have been overcome by including both an iterative and a rearrangement 
procedure that can be requested by the user. However, these options are not provided in this version of 
SAGE even when TWOD is specified, and interested users should contact the authors to obtain a copy of 
SAGE2D and its documentation. For 3-D grid structures, maintaining periodicity at matching planes is a 
more difficult task. Since the initial grid and flow field on matching first and last planes are identical, setting 
the plane torsion parameter (AT) to zero produces matching adapted grids on these planes. In this case, the 
adaption of intermediary planes has no influence on subsequent planes. 

1.11 Appendices 


1.11.1 Appendix I: Derivation of A and B. A and B provide self-adaptiveness to the adaption 
scheme. These two parameters ensure that all repositioned grid nodes maintain grid spacing to within the 
user-requested mesh-size limits (A s MIN and A s MAX ). 

1.11.1.1 Calculation of A. We wish to relate A to the input values of ks MIN and As mx . 
A is constant throughout the entire mesh, and hence the 1-D relationship given in Eq. (2) holds. From the 
original definition of f in Eq. (9), f max = I and f min = 0 ; therefore, from Eq. (3) we have iO max = \ + A and 
(O mm -l . From Eq. (2) (rewritten as As, = K/ to,) we can see that the minimum As will occur at K / (O max . 
Similarly, the maximum As occurs at K/ co m ,„. We therefore wish to set A s MIN = K /( I + A) and A s MAX = K . 
These can be solved simultaneously: by eliminating K we get the expression given in Eq. (10) 

^ MAX _ | 

^ MIN 


1.11.1.2 Calculation of B. This calculation is more complex. B is found by an iterative 
procedure and will change for each j line. The objective is to determine which value of B will give the 


minimum computed As, equal to the requested minimum, 
minimum As,. (An example of a plot of B vs. Av m ,„ is 
shown in Fig. 12). We need to find B at the requested 
As MIN . To do this, we assume an initial value of 
B(=\. 0), evaluate to, and then solve for the new As 
using Eq. (5). Although Eq. (5) is true only for the 
initial line, it is assumed here to hold for all j in order 


As 


MIN' 


For each value of B there exists a 


A s MlN - minks 


(n)\ 


IS 


to simplify the B calculation. If 

small, an acceptable value of B has been found. If not, 
a new B is computed from B (n+]) = B (n) + AB {n \ co is 
reevaluated, and the procedure repeated. 

B (n) can be found from the definition 

■ min As ( t n) ) 


-^-min(As )= lim ^ SsilN — 

dB 3/^o AB (n} 


( 21 ) 


As mentioned in the calculation of A , we know 
that As t is a minimum when co. = l-f-,4, and by 
substituting this in Eq. (5) and differentiating, we 
obtain 
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— mln( As: ) = ■?**■ , 
dB ' 1 + AdB 


a r 1 ^ 


W 


( 22 ) 


where 


We know that 


9 = S 


/=! ®/ 


a ( I 


i a<p 


a^<pj tp 2 as 

Hence, the next step is to evaluate dty/dB. With this definition of <p, we can take the summation sign out 
of the differential and obtain 

a 


dB 


It d 

y— =y— 

^coj £[dB 


\,/=i j 


V®// 


=i- 


1 + Af, t 


n. rB 






/=i 


Equation (22) can now be solved, and after substituting for 

«. ( , 


/=i v®/y 




from Eq. (5), we obtain 


• / a \ >4(1 + A) r . iE'V’ f( B ^ 0 Sfi 

min( Ay, ) = [min( Av, )] \ 

Vr /-l ®/ 


as 


Finally A B can be obtained from Eq. (21), giving B . 

1.11.2 Appendix II: The intersection of a 3-D vector, v, with the arc-length vector, s . This 
technique is used to obtain the b vectors used in the orthogonality terms and to find the location of s' and 
s* . For example, to determine s' we need the segment in which the torsion vector t from the grid point 

(ij-l,k) intersects the j line. To evaluate the direction cosines of the b vector, we need to find the segment 
/— >/ + l along they line that contains the intersection of a vector acting from (ij-l,k) that is normal to 
the j segment. For s* and b* , we are concerned with the line j in the plane k-1 . In all cases, the technique 
for finding / is the same. 

Figure 5 shows the case of the torsion vector t , computed for the point (ij,k) and emanating from 
the point ( ij-l,k ) at D. This vector t most likely will not lie in the plane described by ADC (or the 
equivalent / plane) and thus its projection onto the plane is required. The value of s' in the solution 
equation lies at A\ the intersection of the projection of i onto the ADC plane and the arc-length vector, 
s . The vector n (=n K i +nj + nji) is the unit normal to the plane ADC and is computed by the technique 
described in Appendix III. Both [AA'jand |DA'| are required in the code. 

From vector addition, we can write 

\AA'\s = AD+\DN\i - |M4'|/j (23) 
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) ,A A, 

Since the coordinates of A and D are known, AD(=aJ +a y j + a t k) can be evaluated, hence 
|AA'|. |D/V| and |/VA'| are the three unknowns. By equating coefficients of i, j, and k in the vector 
Eq. (23), we can obtain a set of three equations with three unknowns: 

|AA'|j r = a x + |DAf|/ t - \NA'\n x , \AA'\s y = a y + \DN\t y - \NA'\n y and |AA'|j z = a, + |£W|/. - ]/VA'|n. 

These equations can be solved for |AA'| by computing the relevant determinants. 

\DA'\ is also found by utilizing vector addition. We have 

DA' (= \DA'\ DA') = \AA'\S - AD 

and, since the lengths of the left- and right-hand sides are equal, 

\DAf =(\AA'\s x -a x ) 2 +(\AA'\s y -a y ) 2 +(\AA'\s z -a z ) 2 

1.11.3 Appendix III: Definition of normal vectors. The vector normal to a plane is required 
for the calculation of vectors u, u, b, and b* . This appendix describes the derivation of a general unit- 
vector normal, n, at a given point (ij,k) and normal to the plane ik passing through the constant j—J line. 
This specific orientation is analogous to u. Figure 8 shows such a vector acting at D that represents the 
normal to the shaded surface. In the figure it can be seen that there are four possible planes passing through 
(i,j-l,k)\ GDF, FDM, MDQ, and QDG. Any vector normal required in the code is calculated as the average 
of the normals to each of these four planes. (When the node (ij-l,k) is on or close to a boundary, only one 
or two of these planes will exist.) 

The normal to a plane is derived from the cross product of two vectors lying within the plane, w ith 
care taken to ensure the correct order of the operation (this is why the input grid should be organized as a 

right-handed system). In the example shown in Fig. 8, the four cross products are DFx DM, DFx GD, 

— > — > — > — > 

QDxGD, and QDxDM . The required unit normal is thus the unit vector representing the sum of these 
four vectors. 

This picture gives an idealized view of the defined planes. In reality, some of the required lines have 
already been adapted and some have not. In Fig. 8, the point F could be significantly different from D, 
since D is an already-adapted point and F is still part of the initial grid. The code solves this problem by 
forming a block of data around the current j line that includes all / for lines j — \ ^> j + \ for planes 
A: — 1 => A: + 1 . Lines that have not yet been adapted (e.g .,/+/ at k or j-1 at k+1) are proportioned with 
respect to the s t at j; this relocates the points for a smoother computational effect. 


16 



2 . SAGE USER GUIDE 


2.1 Overview 


The SAGE code is based on the self-adaptive grid method developed by Nakahashi and Deiwert 
(1985). This guide contains information that will enable the user to run the code with little knowledge of 
the mathematical concepts employed in the development of the adaption technique. Included is a detailed 
description of the input-control parameters, along with routine descriptions and nomenclature. The code 
stands alone and has been run on many computer systems. Users wishing to understand and/or amend the 
code can find the details of the mathematical background in the first section of this document. 

The first step in the code is to read two data files: one that contains the coordinates of the grid 
(x,y,z) and another that contains the corresponding flow-field variables, q. It is assumed that both these files 
are in the format (given on the next page) associated with the plotting software package PLOT3D. A 
single-grid, finite-difference, 3-D format is assumed unless re-specified by the user. The next step is to 
adapt the grid with respect to the user-supplied input-control parameters. Adaption takes place as a 
sequence of one-directional adaptions, with the input-control parameters determining the order of 
adaption. Each plane (i.e., 3-D surface) is adapted in one direction only before stepping to the next plane. 
Figure I (see page 2) shows a small section of a 3-D grid: adaption in the i direction (stepping in j) has 
already been performed on the first plane. The code then steps to the next plane and performs the same 
directional adaption on this plane. A 3-D adaption is performed by a sequence of adaptions: once the grid 
has been adapted in one direction, the adaption can be repeated with a new parameter-input set describing 
another direction. (It should be noted that different orders of the adaption direction will produce different 
adapted grids.) When all adaptions are complete, the code sends the redistributed grid points and the 
corresponding interpolated flow variables to two new data files, also in PLOT3D format. 

The analysis that generates the algorithms used in the code is given in detail in the first section of 
this report. Briefly, the redistribution of points is controlled by parameters related to torsion and tension 
springs, and by the maximum and minimum allowable grid spacings. Along each grid line, (x,y,z) is 
transformed to a I-D arc-length variable, s. If we assume that adaption (within a plane, k) steps in the j 
direction, the tension spring constants, co, are evaluated at each point ( i,j,k)\ i.e., CO- =f(s t - f k ) and are a 


function of the gradients of the user-chosen flow-field variables. The torsion parameters (shown in Fig. 2) 

Sj : k _ 2 ) are the link between the current line and the 


Si 




't, = f( hj.k ■ hj-%k ) and V, = fihj.k . 

previously adapted lines. These are the parameters that maintain the integrity of the grid by controlling 
straightness and orthogonality within the plane and between planes. With the evaluation of these variables, 
a system of (n t -2) equations (given in Eq. 8) is developed for the current j line and solved as a tridiagonal 
system. Once the adaption is complete for the current j line, the code steps to the next j line (either 
forward or backward) and repeats the same procedure. At the end of the current k plane, the code moves to 
the next plane and the process is repeated. 


There are eight major steps taken in the code: 

1. Input of three data files: initial grid, flow-field solution, and user-control parameters 

2. Initialization and reorganization of data for computational purposes 

3. Adaption along the start line on the start plane with the 1-D technique 

4. For each subsequent j line on the initial plane: 

a. computation of variables that define the torsion and tension springs (co and x), 
and hence coefficients of the tridiagonal matrix defined by Eq. (8) 

b. iteration to find new values of s, and hence (x,y,z) 

5. Repetition of step 4 until all lines are complete on this plane 

6. For each subsequent k plane, computation of the torsion spring \\f for inclusion in the 

coefficients of the solution Eq. (8) 

7. Repetition of steps 4 through 6 until all planes have been adapted 

8. Output of new grid and interpolated flow- field files in the original format 
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2.2 Execution of the SAGE Code 


The following is an example of the command file to run the SAGE code on a UNIX system. SAGE 
is written in FORTRAN and is self-contained; sage is the executable module. 


cp xyz.gr d fort. 7 
cp q.fun fort. 8 
sage < sage.inp 
cp fort. 10 xyz.out 
cp fort. 11 q. oat 


! copy grid file to unit 7 
! copy solution file to unit 8 

! run SAGE code with sage.inp containing user control parameters 
! name the output adapted grid file 
! name the interpolated function file 


where sage.inp is in namelist format (SNAMEL). {Note: Some systems may need the READ(5,NAMEL,...) 
statement in subroutine INITIAL amended to READ(5,NML=NAMEL,...)}. 


The remaining files are in PLOT3D formatted (ASCII) or unformatted format: 

xyz.grd contains the initial grid points (stored as a right-handed coordinate system) 
q.fun contains the flow- field variables to which the grid is to be adapted 
xyz.out contains the adapted grid points 

q. out contains the flow- field variables interpolated on the adapted grid 

In addition, an output message file is associated with unit 6, which usually defaults to the user’s output 
device. For other operating systems, the user must appropriately assign the six input/output files. 

The following are the read statements for single-grid 3-D PLOT3D unformatted input files: 
xyz.grd: 

READ(7) 

READ(7) 


q.fun: 

READ(8) 

READ(8) 

READ(8) 

where 

IMAX=number of points in the / direction of the grid file 
JMAX=number of points in the j direction 
KMAX=number of points in the k direction 
NDIM=number of Q variables (default=5) 


TMAX,JMAX,KMAX 

(((X(I,J,K),I= I ,IM AX),J= 1 ,JMAX),K= I ,KMAX), 

(((Y(I,J,K),I= 1 ,IMAX),J= 1 ,JMAX),K= I ,KMAX), 

(((Z(I,J,K),I= 1 ,IMAX),J= 1 ,JMAX),K= I ,KMAX) 

IMAX,JMAX,KMAX 

FSMACH,ALP,RE,TIME 

((((Q(I,J,K,N),I=I,TMAX),J=1,JMAX),K=1,KMAX),N=1,NDIM) 


Note: The read statement for formatted (ASCII) files is the same except READ(7) is replaced by 
READ(7,*) etc. 

As seen in the above format, five flow-field variables are expected in the Q File. PLOT3D 
preassigns p, pa, pv, pw, and e, but since the SAGE code requests only the index of the function, any 
variables may be stored. Note that it is possible to handle any number of flow-field variables by changing 
the value of NDIM in the parameter statement at the beginning of each subroutine and recompiling. Also 
contained in the parameter statement are the grid dimensional variables ID, JD, KD, and IMX. These are 
used to define the required maximum dimension of the grid arrays. Because of the internal switching of data, 
these values may not coincide with IMAX, JMAX, and KMAX. If the assigned code dimensions are too 
small, a message will be sent to the user stating the minimum dimension requirements for the current 
application. If the grid is in multiple grid format, an additional record occurs at the beginning of each file 
stating the number of grids, and the grid dimensions are supplied in dimensioned arrays. FSMACH, ALP, 
RE, and TIME are not used in the code and may contain dummy values. They are part of the PLOT3D 
package and are displayed on the output plots. 


2.3 User Input Parameters 

The file sage.inp is the user-supplied, input-parameter control file. The grid adaption is based on 
the user’s choice of these input parameters, which are listed and briefly described below. This is followed by 
a more complete explanation of each parameter. The input file sage.inp uses the namelist format, since 
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generally only a few of the input parameters need to be changed from the default value set by the code. 
These default values are shown in parentheses in the list given in section 2.3.1. If more than one adaption 
pass is to be made (for example, a two-directional adaption or multiple grids), the subsequent adaptions can 
be made on the adapted grid by linking up to 10 sets of namelist inputs within the same sage.inp file. For 
multiple grids, these sets can also contain multiple export/import passes. 

Before describing the input parameters, some terminology needs to be clarified: 

The term physical domain is used to reference the complete grid defined by the input grid file, i.e., 
the grid bounded by IMAX, JMAX, and KMAX. The adaption domain is the part of the grid, as defined by 
the input-control file, that is to be adapted, i.e., (1ST, TEND), (JST,JEND), (KST,KEND). These two 
domains can be equivalent. The direction /, y, or k refers to the direction of the grid coordinates as defined 
by the order in which they are stored in the grid file. The first index (normally containing x) is named the / 
direction, the second index is the j direction, and the third index the k direction. This implies that if data 
happened to be stored as {z,x,y) instead of (x,y,z) 9 / would represent the z direction. Regardless of the order, 
the data must be stored as a right-handed coordinate system. The adaption direction is used to define the 
1-D line along which adaption (redistribution of points) takes place. In the analysis and in the descriptions 
below, this direction is always / for convenience, but the user may request i, j , or k . There are two stepping 
directions : one within the plane, defined as j in the analysis, and the other in the direction of plane 
marching, defined as k. Although the analysis and descriptions in this report assume this particular order of 
adaption and stepping, the code makes no such a priori assumptions, since the order is controlled by the 
user’s input-parameter file. Reference is also made to grid boundaries . Side-edge boundaries refer to the start 
and end points at the edge of the adaption line (assumed to be / in the analysis). These are the edges that are 
controlled by the NEDGE parameter. Marching boundaries are the start and end lines of the stepping 
(within the plane) direction (assumed to be j). Adaption close to these lines are affected by the parameters 
ORTEIS(l), ORTHE(l), and MGSTEPS, as well as by some internal controls. Planes juxtaposed to start and 
end planes are affected by ORTEIS(2), ORTHE(2), and MGPLS. 

2.3.1 Parameter control file, sage.inp . The input parameters, with their default values (in 

parentheses) and short descriptions, are described below. They are not in alphabetical order, but grouped by 
category or affinity. Reminder: $NAMEL must begin in column 2. 


$NAMEL 



FORMI 

(FALSE) 

set to .true, if input PLOT3D files are formatted (ASCII) 

FORMO 

(FALSE) 

set to .true, if formatted output PLOT3D files requested 

1ST 

(I) 

first adaption line in / direction 

IEND 

(IMAX) 

last adaption line in / direction 

JST 

(1) 

first adaption line in j direction 

JEND 

(JMAX) 

last adaption line in j direction 

KST 

(1) 

first adaption line in k direction 

KEND 

(KMAX) 

last adaption line in k direction 

IJPLANE 

(TRUE) 

the adaption surface lies in the (ij) plane, with 
plane stepping occurring in the £ direction 

JKPLANE 

(FALSE) 

the adaption surface lies in the (j } k) plane, with plane stepping 
occurring in the / direction 

IKPLANE 

(FALSE) 

the adaption surface lies in the (i,k) plane, with 
plane stepping occurring in they direction 

ISTEP 

(FALSE) 

=.true. for stepping in the / direction within the plane: 
e.g., it cannot=.true. if JKPLANE =.true. 

JSTEP 

(TRUE) 

-.true, for stepping in they direction within the plane 

KSTEP 

(FALSE) 

=.true. for stepping in the k direction within the plane 

INDQ 

(1) 

index of adaption flow-field variable, q\ default of 1 => p 

IQ(8) 

(0) 

enables a combination of q variables to drive the adaption 
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RDSMAX 

(2.0) 

RDSMTN 

(.5) 

CLAM(I) 

(.01) 

CLAM(2) 

(.0001) 

CT(!) 

(•5) 

CT(2) 

(.5) 

NEDGE 

(0) 

MG1 

(0 or 4) 

MG2 

(0 or 4) 

INTER 

(2) 

NFILT 

(2) 

MGSTEPS 

(0) 

MGPLS 

(0) 

MARCH 

(FALSE) 

MARCHPL 

(FALSE) 

ADD 

(0) 

LSTADD 

(1ST) 

LENDADD 

(IEND) 

SUB 

(0) 

LSTSUB 

(1ST) 

LENDSUB 

(TEND) 

REMOVE 

(0) 

NOUP 

(FALSE) 

SAVE 

(TRUE) 

ORTHS(2) 

(TRUE) 

ORTHE(2) 

(TRUE) 

GEOM 

(FALSE) 

QFUN 

(TRUE) 

NOQ 

(FALSE) 

LNSING 

(0) 

PLSING 

(0) 

TWOD 

(FALSE) 

FV 

(FALSE) 

MGRID 

(0) 

EXPORT 

(FALSE) 

IMPORT 

(FALSE) 

MPLANE 

(1) 

IS,IE 

(0,0) 

JS,JE 

(0,0) 

KS,KE 

(0,0) 

MVBOUND 

(0) 

RECLUST 

(0) 

NSM 

(10) 

NSM1 

(10) 

NOSESM 

(0) 

DSW 

(0) 

DSE 

(5.0) 


relative maximum allowed As:As max > 1.0 
relative minimum allowed As:As MfN <1.0 

=X, coefficient of torsion parameter T, 1 <X< lO^ 6 
= AT, coefficient of plane torsion parameter, \jr; same range as X 
proportion of "straightness" to "orthogonal" for torsion vector T 
as CT(1), but proportions plane torsion vector, /* 

override control on side-edge adaption: 
l=both edges, 2=start edge, 3=end edge 

used w ith NEDGE: number of points to merge start-side spacing 
used with NEDGE: number of points to merge end-side spacing 
order of interpolation: 2, 3, or 4 
number of passes to filter (smooth) q and co 

number of merging lines for within-plane torsion controls 
number of merging planes for between-plane torsion controls 
=.true. to extrapolate end adaption line throughout remaining lines 
=.true. to extrapolate last adapted plane throughout remaining planes 

= 7 ? to add n points between each node in selected range 
lower limit of range for adding points (i.e., when ADD>0) 
upper limit (if ADD=0, values are ignored) 

—n to delete n points between each node in selected range 
lower limit of range for point deletion (i.e., SUB>0) 
upper limit (if SUB=0, values are ignored) 
removes requested number of points from outer grid region 

= true. if no adaption required (e.g., for moving outer boundary ) 

=. false, to suppress output of data files 

=. false, to remove orthogonal constraint at start boundaries 
=. false, to remove orthogonal constraint at end boundaries 

=.true. to include geometry at wall boundaries in adaption variable 
=. false, to remove f(q) from adaption variable 
(used with GEOM=.true. only) 

=.true. if the q file is not used 

-n if the line n is common to all adaption planes 

=n if adaption plane n is collapsed to a line 

=.true. if datasets are 2-D (see section 2.4.1) 

=.true. if q file in finite-volume format (2.4.2) 

adaption grid number if multiple grids 

denotes an export plane transfer, not an adaption pass, multigrid only 

denotes an import plane transfer, multigrid only 

defines plane number, used for muftipfe grid export and import only 

defines i direction domain for plane transfers 

defines j direction domain for plane transfers 

defines k direction domain for plane transfers 

nonzero to move the outer boundary 

method to redistribute points, with or without boundary move 

smoothing filter in the line-stepping direction for outer boundary 

smoothing filter in the plane-stepping direction for outer boundary 

smoothing across planes with singular line at nose 

wall spacing, only for RECLUST>0 

outer edge spacing, only for RECLUST=1 or unax 
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DSW1 (0) 

DSW2 (0) 

DSN (5) 


only for DSW=999, wall spacing at j=l 

only for DSW=999, wall spacing at j=jmax 

percentage to expand when matching outer boundary to shock 


First adaption Although many parameters have just been described, generally only a few are used 
for each adaption. The best technique, even for experienced users, is to retain as many of the default 
parameters as possible, view the results and then adjust some parameters if necessary. For some 3-D 
problems, it is clear which plane should be the stepping plane so the PLANE parameter, along with the 
STEP parameter, can be initially chosen. Note that only one PLANE parameter and one STEP parameter 
are needed for one adaption pass. Also, a flow-field variable should be chosen for TNDQ that shows the flow 
features most clearly . To find appropriate with in-plane parameters, adapt on only one plane and/or turn 
off the connectivity between planes (i.e., CLAM(2)=0) to see the set of 2-D adaption planes. If there are 
any lines common to all planes or if a plane is collapsed to a line, PLSTNG and/or LNSING must be set. If 
the first adaption pass is hot acceptable, try increasing the ratio between RDSMAX and RDSMIN, 
decreasing CLAM(l), and/or setting NEDGE=1. Since many of the parameters have interdependent effects, 
it is better to change only one or two of the parameters at a time. 

2,3.2 Explanation of user-supplied input parameters. The following is a detailed explanation 
of each of the input parameters; they are listed in alphabetical order. 

ADD When this is nonzero, points are added between adjacent mesh points within the requested 
range (see LSTADD, LENDADD). For example, if ADD=2, two points will be added between each 
consecutive grid point. Adding occurs only in the adaption direction and not in either of the stepping 
directions. Ensure that the added points do not cause the coded array dimensions to be exceeded. ADD and 
SUB can be used in the same pass to move points: note that the ADD will occur before the SUB, so the 
range on the SUB parameter needs care. 

CLAM CLAM(l)=?l and CLAM(2)=A* define the magnitude of the torsion parameters t and \| /, 
respectively. As the values of these parameters decrease, more points will be pulled into the high-gradient 
regions, at the possible expense of grid smoothness. 

CLAM(I) controls x, the torsion parameter within the plane; its order of magnitude can lie between 10 -6 
and 1. A value of zero produces a set of independently adapted lines, possibly generating crossed grid lines. 
As A increases, the grid becomes smoother but less adapted. 

CLAM(2) controls the magnitude of \j /, the torsion between planes. It has the same range of values as 
CLAM(l); however, if it is zero, the adapted grid may still be acceptable. For periodic planes (i.e., if the 
first and last planes are the same), setting CLAM(2)=0 is necessary to prevent a discontinuous grid at the 
juncture. 

CT CT(1)=C, and CT(2)=C*; they represent the direction of the torsion vectors 7 and 7 * 
(whereas X and X* are their magnitude in these directions). They have the range of 0 < C / 5 C * < 1.0, where a 
value of zero emphasizes orthogonality and a value of one emphasizes straightness. The default of .5 places 
the torsion vectors halfway between. This value is suitable in most cases, but it may cause problems when 
side boundaries are concave or w hen adapting on already adapted planes. 

DSE This is used only when RECLUST=/7«ax (note: RECLUST=1 assumes RECLUST=/>/?ax) and 
provides the outer edge spacing. The Vinokur algorithm needs the first and last edge spacing along the 
recluster line. DSE is the multiple of the average mesh size As, and is defaulted to 5. Hence the last-edge 
mesh size input to the algorithm is DSE *s mtx /(npts-\) where s max is the new s nwx if the outer boundary 
location has changed. If 3>RECLUST</w<zr, the code computes the outer edge spacing to match across the 
boundary with the existing grid spacing and DSE is not used. 

DSN When MVBOUND is based on a shock distance or a flow contour, the code finds the 
requested flow feature, s slufi , k , and places the new boundary point outside'of it. DSN (defan It=5) determines 
the new boundary location as s shock + DSN*s, mjv /I00. , where s mflLT is the old s max . 
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DSW This is the first mesh spacing at the inner boundary' for the Vinokur algorithm (RECLUST* 
0). If no value is input, the original As, is used. DSW can be very small if necessary', but if too small, there 
may be insufficient latitude to permit the algorithm to reach an appropriate solution downstream, creating 
a jagged grid. This problem may be resolved by providing a first and last wall spacing (DSW1 and DSW2). In 
this case, set DSW=999, and input DSWI and DSW2. The code will compute DSW for each line. 

DSW1, DSW2 These two parameters replace the single DSW and may resolve the jagged problem 
mentioned above. A very' small value of DSW may be appropriate at the first active line (j=l) where s max is 
small, but if the grid expands rapidly, the value of s mux downstream is too large to accommodate the same 
small wall spacing. In this case, DSWI will be used at j—1, and DSW2 will be used at j—jmax and the 
intermediate values will be interpolated. 

EXPORT For multiple grids only. When EXPORT is set to .true., the user-supplied parameters for 
this adaption pass describe the current location of a plane of data to be transferred from one location to 
another, most probably in another grid. This parameter set must be followed by an IMPORT parameter set 
that describes the plane’s destination. 

FORMI If equal to .true., the PLOT3D grid and Q input files (assigned to units 7 and 8) are 
assumed to be ASCII (formatted). Default is unformatted. Note that if SAGE is compiled as r8 (double 
precision), then only formatted files can be used. 

FORMO If equal to .true., output files will be formatted (the unit numbers depend on the number 
of adaptions and multiple grid rereads). Default is unformatted. 

FV This is set to .true, if the q file is in finite-volume format (i.e., the q variables are evaluated at 
the cell centers and not at the grid points). See section 2.4.2 for a more detailed description of the finite- 
volume option. 

GEOM This parameter should be used when a wall boundary is defined by high surface gradients 
and the standard grid redistribution has moved points in such a way that the original shape has been 
deformed. When GEOM is set to .true., the code will add the surface curvature function to the flow-field 
gradient function in the wall boundary regions. Points will thus be maintained in or redistributed into regions 
of high surface curvature as well as into regions of high flow-field gradients. The contribution of the 
geometry function will proportionally decrease away from both boundaries, with the internal lines 
controlled by the flow field only. If GEOM=. true, and QFUN=.false., adaption will be to geometry gradients 
only, for all grid lines. 

IJPLANE, rKPLANE, JKPLANE These parameters define the plane on which adaption takes 
place. By default, they also imply the direction of the stepping plane. The input parameter file needs only 
one of these parameters to be set to .true., and the code will automatically assign .false, to the other two. 
IJPLANE=.true. indicates that the plane represented by (i,j) will be the adaption plane and that plane 
stepping will occur in the k direction. Whether k is a forward or backward step will depend on KST and 
KEND (i.e., if KST > KEND then backward stepping will occur). Similarly, IKPLANE=.true. indicates that 
the adaption plane contains the (i,k) directions and that j is the plane-stepping direction. Finally, 
JKPLANE=.true. refers to the plane containing the points (j,k) and i is the plane-stepping direction. 

IMPORT For multiple grids only. An input parameter list containing IMPORT=.true. must 

immediately follow a parameter set containing EXPORT=.true. The IMPORT list describes the destination 
(i.e., the receiving plane) of the plane of data defined in the EXPORT=.true. input-parameter set. The use 
of MGPLS is important for the subsequent adaption of the import grid. 

INDQ This parameter indicates which of eight possible flow-field variable(s) will drive the 
redistribution of grid points. From the standard PLOT3D format, five options are available: 

1 — » density p 

2 — > x-momentum, pu 

3 — > y-momentum, pv 

4 — > z-momentum, pw 

5 -» stagnation energy, e 

Three more options are available by setting FNDQ=6, 7, or 8 
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6 -> pressure, p 
1 Mach number, M 
8 — > temperature ratio, T 

Pressure, Mach number, and temperature ratio are computed using the ideal gas relationship and 
assumes the Q file contains the standard variables. (The code actually assigns pressure to NDIM+1, Mach 
number to NDIM+2, and temperature to NDIM+3, so if the user has changed the value of NDIM to 
accommodate extra flow-field variables, INDQ must reflect this change.) Note that INDQ=4 is not 
available for 2-D datasets. The user wilf normally choose to adapt to the flow-field variable that most 
represents the flow features. However, if different variables demonstrate different features, it may be 
advantageous to combine them to bring out all the features on the adapted grid. In this case, set INDQ=0 
and input values for IQ. 

INTER indicates whether to use a linear (INTERS), a quadratic Lagrange polynomial (INTER=3), 
or a cubic spline (INTER=4) scheme for interpolations. Interpolation is used throughout the code; for 
example, the q values in the output function file are interpolated at the new adapted grid points. Linear 
interpolation will usually provide the appropriate result. 

IQ is an array of eight (or NDIM+3) integer values that are used only when INDQ=0. They allow 
the user to modify the adaption variable to a combination of variables. The order of IQ is consistent with 
the order of the flow-field variables in Q. The value of an index is the proportion that the corresponding 
variable will contribute to the final adaption variable. For example, TQ( 1 )= 1 ,IQ(7)= 1 [i.e., 

(1,0, 0,0,0, 0,1,0)] will produce an adaptive function of — and TQ(I)=1, TQ(3)=1, IQ(6)=3 will 

2 ds ds 

produce 1^ + 1^ + -^. Obviously, IQ=( 1 ,0,0, 0,0, 0,0,0) is the same as INDQ=1. 

5 ds 5 3s 5 ds 

1ST, IEND contain the indices defining the first and last boundary lines of the adaptive domain in 
the / direction. Similarly, JST, JEND define the domain in the j direction and KST, KENT) define the 
domain in the k direction. These variables define the limits of the adaption domain and must lie within the 
input grid boundaries defined for the physical domain, i.e., 

1 < 1ST, TEND < IMAX; l< JST, JEND < JMAX; and 1 < KST,KEND < KM AX 
Forward and backward stepping are also controlled by these parameters. If plane stepping is in the k 
direction, setting KST > KEND will produce backward stepping. Similarly, if stepping within the plane is in 
the j direction, setting JST > JEND will produce backward stepping. The reversing of the data is handled 
internally and is imperceptible to the user. Reversing either of the stepping directions will completely 
change the resulting grid, since the redistribution along the initial line and plane will be different, as will the 
connecting torsion springs. Reversing the order of the adaption direction (i.e., / in the default case) should 
have no effect on the solution, since the solution along a line is independent of the order of points. 
However, it can be used to redefine start and end points on a line if necessary, and for meshes that contain 
very large and very small A s /9 numerical accuracy may be influenced by this adaption direction. 

ISTEP, JSTEP, KSTEP These are used in conjunction with the PLANE parameters described 
above and define the marching and adaption directions within the defined plane. Only one of these three 
parameters should be input and set to .true., and the code will assign .false, to the other two. If 
IJPLANE=.true., then only ISTEP or JSTEP can be true (KSTEP must be false). If ISTEP=.true., stepping 
occurs in the / direction and thus the adaption direction will be j. Points will be adapted along each constant 
/ line and stepping will occur to the next / line, forward or backward, depending on 1ST and TEND. If 
JKPLANE=.true., only JSTEP or KSTEP can be true and similarly, for TKPLANE=.true., only ISTEP or 
KSTEP^.true. will have any meaning. The code puts out an error message if these inputs are inconsistent. 

IS, IE, JS, JE, KS, KE These parameters are used only for an export or import transfer process in 
a multiple-grid file. These variables define the range of the transfer domain within the plane, MPLANE. 
The PLANE parameter (e.g., IJPLANE) defines which coordinate plane is being transferred. As an 
example, if IKPLANE= true., then IS, IE, KS, and KE are used to define the domain within the transfer 
plane. If they are omitted, it is assumed that the entire plane is being transferred. 

LNSING In some 3-D grid types, planes emanate from a common grid line. If the chosen set of 
adaption planes includes this common line, then this line should be adapted only on the first plane and not 
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on subsequent planes. This adapted line is then placed in the first line of all planes. Although the input 
option is LNSING=n, where n is the line number, it is likely that n = l . 

LST ADD, LEND ADD These are input only if ADDzD and if only a portion of the grid is to be 
expanded. If they are not input and ADD^O, the entire adaption domain (not physical domain) is assumed. 
If ADD=0, their values are ignored. 

LSTSUB, LENDSUB These parameters are input only if SUB^O; they define the limits in which 
points are to be removed. If they are not input, then the entire adaption domain is assumed. 

MARCH This parameter refers to stepping within the plane. If the last line to be adapted ( j end ) is 
within the physical grid boundary (i.e., j end < j„ utx ), a sharp discontinuity will occur between the last adapted 
line, j end , and the nonadapted line, j end+i . Setting MARCH=.true. causes the remaining lines within the 
plane (i.e., j end+ , j mai ) to be realigned so that they are proportional to the last adapted line. This 
realignment will be performed for every plane. 

MARCHPL This parameter refers to the plane-stepping direction, and can be used independently 
or in conjunction with MARCH. If the last adaption plane is within the physical boundary (i.e., k end <k max ), 
each line in each subsequent plane will be proportioned with respect to the adapted lines in the k end plane. 

MG1, MG2 When NEDGE is nonzero, an override mesh spacing is computed at the requested 
boundaries (either first, last, or both). This edge-point spacing, which is not a function of the adaption but 
of the initial grid, is merged into the three adjacent points to produce a smooth transition. The default 
value of MG1 (for /„) and MG2 (for i a , d ) is zero when NEDGE=0, but four when NEDGE is requested. The 
user has the option of overriding this value, setting it to any other integer value. This can be used when 
adapting a boundary layer; when either MG is increased, the dense edge spacing is maintained over a larger 
region. For example, MG 1 = 10 would merge the edge spacing into the next nine cells. 

MGPLS This input variable is analogous to MGSTEPS described below, but applies to planes. If 
the first adaption plane (k sl ) is internal to the physical domain, MGPLS=« permits the user to gradually 
bring in the effect of the adaption parameters to produce a smooth transition across the start plane. No 
adaption will take place on the first plane (&,.,), and after n planes, full adaption occurs with the adaption 

parameters C* , X, etc., equaling their input value. This is an especially important feature for grids that 
have an initial distribution on the wall boundary that defines a physical shape and cannot be changed. It is 
also useful for retaining matching multiple-grid boundaries (see IMPORT). 

MGRID This parameter is only used for multiple-grid files and indicates which grid to adapt. Do 
not use MGRID=1 for a single grid: setting MGRID to nonzero automatically implies a multiple-grid file. 

MGSTEPS (m ) provides continuity when the first adaption line on each plane is internal to the 

physical boundary. Inputting MGSTEPS =/7 tells the code to start with no adaption on the initial line (i.e., 
retain the original distribution on the j sl line) and to linearly increase the adaption effect until, after n 

lines, full adaption occurs. At this point, C„ X, etc., will coincide with their input values. If MGSTEPS=l, 
no adaption will be performed on the first line, but full adaption will occur on the second line. Figure 9 
show r s an example with m = 5, and section 1.4.1 describes the parameter in detail. 

MVBOUND This is the input parameter that invokes the outer boundary movement. There are 
three options: 

(1) MVBOUND=±p, where p is a percentage of the input line length. For example, 
MVBOUND=-IO. will create a grid where s max along each line is 10% smaller than the initial grid. 

(2) MVBOUND=999. will move the outer boundary based on the shock location s shllck of the 
specified variable, INDQ. The code will look for a sudden change in gradient of the flow variable and add a 
given distance to this shock location to create a new outer boundary. 

(3) MVBOUND>1000. will allow the boundary to parallel a given contour value found in the Q file 
(defined by INDQ). The contour value is MVBOUND- 1 000. 
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MPLANE Used only for multiple grids during an import or export process. MPLANE is an integer 
defining the transfer plane. The PLANE parameter is used (e.g., IJPLANE=.t.) to indicate the coordinate 
direction of MPLANE. 

NEDGE is a flag that requests an override on the computed side-edge boundary spacing. Side edges 
occur at i s( and i eml and frequently need special handling. If there are no flow gradients near the edge of the 
domain, the standard adaption algorithm will pull points away from the edge. This may not be a 
satisfactory result, as, for example, when the side edge is internal to the physical grid boundary. 
Figure 10(a) shows a side-edge adaption with NEDGE=0. The flow-field gradients are concentrated in the 
center of the grid, and the first adaption point (i=4) has been pulled far from the boundary line at i sf =3. It 
is clear that it is preferable for the adapted side-edge spacing to be continuous with the juxtaposed spacing 
in the external region. Even if the two boundaries coincide (i.e., i st = 1), the user may prefer a different 
spacing than that computed by the adaption algorithm. In either case, NEDGE can be set and the code will 
try to improve the side-edge spacing. The computed mesh-size override is merged into the next four points, 
but this number can be changed by MG1 and MG2. The result of setting NEDGE=1 is shown in Fig. 10(b). 
Depending on the case, both or only one of the side spacings may need improving. NEDGE=1 requests both 
edges, NEDGE=2 requests start edge only, and NEDGE=3 requests end edge only. 

NFILT is a "filtering" variable that defines the number of passes used to smooth the gradient of the 
input q data and the computed tension parameter, co. The default value of two will generally suffice, but if 
the flow-field variables are discontinuous, it may be helpful to increase the value of NFILT. An increase in 
NFILT can also be used to expand or spread out a very sharp flow feature. 

NOUP This variable is used to change the grid (e.g., ADD, SUB, MVBOUND, RECLUST) without 
performing an adaption. NOUP stands for NOUPdate. 

NOQ If no q file is available, SAGE can still be used to smooth the grid; perhaps to equal spacing 
or to the geometry function, or to recluster, with or without moving the outer boundary. If NOQ = .true., 
no file will be read on unit 8, and instead a constant flow field will be generated internally and there will be 
no interpolated solution file. This is quite different from QFUN = .false, where the q variables are not used, 
but are interpolated onto the new r grid and output. 

NOSESM This variable was added to remove a problem caused by a special grid structure that uses 
a singular line at the nose. Since the smoothing parameters (NSM, NSM1) fix the first and last points, s max 
at the nose will remain fixed. NOSESM will change the nose s max to give a smoother appearance. NOSESM 
is the number of adjacent points used to smooth s max over the nose region. Typically the value should be 1 , 
2, or 3. 

NSM, NSM1 NSM and NSM I use the same filtering option described elsew here (see NFILT and 
the FILTER subroutine). Since the outer edge of a shock may not be smooth (due to initial calculations 
and/or a coarse mesh) these two parameters will help smooth the new boundary, and though the default 
value is 10 in both cases, numbers up to 100 may be appropriate. NSM is used to smooth the i-imax for all 
j in the IJPLANE, whereas NSM1 smooths the cross plane {i-imax for all k in the IKPLANE). 
Occasionally this smoothing parameter overemphasizes an irregularity in the outer surface and should be set 
to zero if this occurs (see Fig. 35 in Section 3). 

ORTHS, ORTHE The code assumes that orthogonality to the marching boundaries is desirable. 
This may not be the case, as, for example, in outgoing flow where the shape of the outer boundary is 
arbitrary. Setting ORTHS(l) = .false, will turn off orthogonality from the first to second adaption lines, and 
ORTHE(l) = .false, performs the same function when approaching the end adaption line. ORTHS(2) and 
ORTHE(2) will similarly affect the plane boundaries. 

PLSING This is a parameter unique to 3-D grid adaption. It stands for plane singularity, and 
implies that a plane n is actually a single line, but is stored as a set of identical lines. The code will not be 
able to compute normals and will certainly "blow up" unless informed of this condition. It is only relevant 
when the adaption plane direction coincides with this collapsed plane. Do not use it when adapting in other 
plane directions. 

QFUN This parameter permits the user to adapt to the geometry function only. When QFUN = 
.false, and GEOM = .true, are input, the coefficient of the flow-field gradient (C <y ) is set to zero for all grid 
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lines (see section 1.5). In this case, the geometry function is computed throughout the grid and drives the 
adaption for all grid lines. 

RDSMAX, RDSMIN control the density of the redistributed points and are the maximum and 
minimum allowable grid spacings. They are input as proportioned values and are changed to physical 
variables internally, i.e., As max X5„ I(K /(/j ; -U and As MIN x s miLt /( //, - 1), where //, is a constant equal to the 
total number of points along the adapted line, and s max is the length of the current adaption line. This 

implies that RDSMAX>1.0 and RDSMIN<1.0. (If both are set to 1.0, a uniform grid will result.) As an 
example, RDSMIN = 0.5 will prevent a converged As from being less than half the average step size. (Since 
many factors influence the distribution of grid points, this control is not absolute.) Note that since the 
adaption along the first line is not influenced by the torsion parameters, this initial line will present a 
clearer picture of the effect of RDSMAX and RDSMIN. 

RECLUST This parameter is used to recluster the grid points, with or without boundary 

movement and is quite different from adaption. If outer boundary movement is requested 
(MVBOUND^O), the default of RECLUST=0 specifies that the new spacing on the line is proportional 
to the original spacing. RECLUST=I will invoke the Vinokur algorithm and provide spacing based on an 
inverse hyperbolic function. If no outer boundary movement is desired (MVBOUND=0), but reclustering 
based on the Vinokur algorithm is required (RECLUST=/i, where 3 > n), then n points will be reclustered. 
This is useful for limiting reclustering to the boundary layer or for entire grid reclustering. Note that this 
grid redistribution option is also useful for removing zero cells (i.e., two grid points at the same location) 
from an initial input grid. 

REMOVE It is possible that an initial grid has unnecessary grid points in the outer region. Once an 
initial solution has been obtained, the user can see that these points are wasted. REMOVE=« will remove n 
points from the end of all adaption lines. NOUP should be set to true if only removing points is required. 
To remove points from the inner region, reverse 1ST and TEND. This will not change the adaption but will 
fool the remove operation. Remember, REMOVE deletes the boundary points, SUB retains the outmost 
one. 

SAVE This is useful when more than one adaption pass is made in the same program run, for 
example, an adaption stepping in they direction followed by one stepping in the / direction. The output grid 
and function files are large, and setting SAVE=. false, on a set of SNAMEL will suppress the output for that 
particular adaption. If SAVE=.true. (default), each subsequent output set of xyz.nut and q.out files will be 
assigned to different unit numbers. As stated in the execution section, the first output set is assigned to units 
10 and 1 1. The second output set will therefore be units 12 and 13, and so on. Note: SAVE=.faIse. cannot 
be used for multiple grids. 

SUB When SUB=/z (*0), points are removed from the adaption line. As an example, if « = /, every 
other point is deleted. If «=2, two consecutive points are deleted between the points retained. Note that the 
number of stepping lines remains constant: points are only removed from the adaption line. To remove 
points from both directions, two passes are required. See LSTSUB, LENDSUB if only a selected range of 
deleted points is to be deleted. 

TWOD If the input grid and function files are stored as 2-D PLOT3D files, this parameter must 
be set to .true. The code will assume that IJPLANE is the adaption plane and JSTEP the stepping direction. 
The user must specify ISTEP=.true. if required. The next section discusses the adaption of 2-D problems. 

2.4 Alternative Grid Types 

2.4.1 Two-dimensional adaption. The SAGE code can accommodate 2-D datasets, and will 
adapt the single plane in the same manner as in the original 2-D SAGE code (Davies and Venkatapathy, 
1989). When the input parameter TWOD is set to .true., the code will read the grid (assuming (x,y)) and 
function files as 2-D files. These datasets will then be internally converted to a 3-D format; the number of 
k planes will be set equal to one, creating a constant r coordinate, and the 4th q function will be shifted to 
index 5. (Note that indices INDQ and IQ retain their 3-D relationship.) Because of this reorganization, no 
special handling of 2-D datasets is required within the body of the code. Datasets are reconverted to 2-D 
form before output. To accommodate larger dimensions, the parameter statement may be changed at the 
beginning of each routine. Since 2-D datasets require only KD=1, ID, JD, and IMX may be significantly 
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increased. However, this change is not made automatically and the user must change the parameter 
statements if necessary. 

2.4.2 Finite-volume grids. The solution file associated with finite-volume applications contains 

q values evaluated at the cell center. Therefore the IMAX, JMAX, and KMAX values on the header record 
are one less than the values given in the grid file. If SAGE finds that the size records disagree and that the 

finite volume option is off, an error message is sent to the user. 

Care must be taken with the boundary (or ghost) cells. SAGE interpolates q onto the internal grid 

points and then sets all boundary values (in the physical domain) equal to the adjacent interior value. After 

adaption, flow values are interpolated back to the cell center. If all planes are not adapted, the q values in 
the final adapted plane will be interpolated as if they are on a 2-D surface since there will be a discontinuity 
between the adapted and non-adapted grid points. The boundary values of the physical domain are again set 
equal to the adjacent values, regardless of whether the entire physical domain has been adapted. It is 
therefore very important for the user to check all the boundary cells in each coordinate direction. 

Related to the handling of the boundary cells, SAGE will only adapt a finite-volume grid if it is a 
single surface (whether defined in 2-D or 3-D) or is a 3-D grid with four or more planes. 

2.4.3 Multiple grids. The multiple-grid format is a single file containing a collection of 
separated grids. Preceding the grids are two header records, one defining the number of grids and the other 
the size of each grid. The associated q file is similarly defined. 

Some complex multiple grids utilize the blanking feature available in PLOT3D. Currently, this 
version of the SAGE code does not handle these blanked regions (see Section 5). Also, if overlapping 
regions are adapted in separate grids, it is the user’s responsibility to interpolate the results. However, 
matching zonal planes can be handled with the plane-transfer feature described below. 

2.4.3. 1 File handling. Since each grid is stored sequentially within a multiple-grid file, 
SAGE reads and copies all grids (and their associated solutions) to the appropriate output files until the 
requested grid (MGRID=/?) has been read. This grid is now adapted (or a plane transferred) and written, 
along with the interpolated flow solution, into the correct sequence in the output files. Finally, any grids 
following the adapted grid are also read and copied to the output files. Subsequent adaption or export/import 
passes will use these output files as input files. If several passes occur in the same computer run, several sets 
of large files could be created and the user is reminded to delete unneeded files. The SAVE parameter cannot 
be used since the code needs to rewind the files for multiple passes. 

An additional set of files will be created if any grid has changed size (e.g., the ADD or SUB option 
has been used). Grid sizes are stored in the header record that has already been written to the output file 
before adaption took place. This header record must be amended to reflect the new grid size. The final 
output files are therefore read in as input files, the header record is corrected, and yet another set of output 
files is created. An output message keeps the user informed of the unit numbers for the final set of output 
files. 

2.4.3. 2 Data transfer using EXPORT and IMPORT. Section 1.7 explains the need for 
data transfer between grids. A set of input parameters is available that controls the transfer of data from 
one domain within a specified plane of a specified grid to a matching domain within another grid. The same 
namelist format is used that controls the normal adaption procedure in SAGE, but two sets are required in 
the user-input file: the first describes the ‘export’ plane and the second describes the ‘import' plane. Here, 
the term ‘export’ means a domain (i.e., a surface, plane, or sub-plane) whose data will be transferred. The 
domain receiving this data is the 'import’ domain. A transfer domain is defined by 

1. the transfer code, either ex port=.t. or import-. L 

2. the grid number, mgrid 

3. the plane direction, (ijphme,ikplcmejkplane) 

4. the plane number, mphme (this will often be the first _or last plane) 

5. the range of points describing the domain (i.e., the start and end points in two directions) 

within the plane 

Notes: (a) For 2-D datasets only one line in a plane is transferred and ijplane =./. and mplane = I 
are defaulted by the code. The transfer domain is described by is, ie, js , je and one of these directions must 
have equal start and end points. If the ‘export’ card is the first input set, remember to include hrod=.t. 
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(b) Multiple transfers can be handled in the same run of the code as well as a combination of 
adaptions and transfers. See the examples section for clarification. 

2.4.4 Blanked grids. SAGEv2B is the only version of SAGE that recognizes the blanking option 
that is found in some PLOT3D grid files as an additional variable, called IBL. Section 5 of this document 
describes the blanking option in detail. If IBL is defined in the PLOT3D input grid file, SAGEv3 can be used 
but the blanking will be ignored and the output grid file will not contain IBL. 

2.5 Output Message File 

The sage. out file is written to unit 6, the normal default for the output screen. It contains messages 
that help explain what has happened during program execution. At the end of each adaption, "ADAPTION 
n COMPLETE" indicates that the program was able to run to completion and that xy.out and q.out files 
have been created. The message "OUTPUT FILES ON UNITS n, AND n 2 " informs the user that the 
output grid and Q files are named fort.n, and fort.n 2 . 

The following are other messages that may be seen (given in alphabetical order) along with a short 
description of their meaning. 

ADD OPTION EXCEEDS DIMENSION (Critical) 

Self-explanatory. Increase array dimensions. 

CANNOT USE INDQ > NUMBER OF Q IN SOLUTION FILE (Critical) 

This message comes from subroutine SHOCKLE. Only variables within the Q file can be used to detect 
a shock, not internally computed values. 

FINITE VOLUME METHOD NEEDS I OR 4+ PLANES (Critical) 

A single plane is treated like a 2-D surface: i.e., 4 points are used for cell-centering. Four planes are 
needed for the finite volume method in 3-D. 

GRID HAS IDENTICAL POINTS AT / AND i+I ON LINE j AND PLANE k: 

USE RECLUST OPTION TO TRY TO REMOVE (Critical) 

Computation of body normals is impossible in cells of zero size. Use the reclust option to remove 
them, and perform the adaption again. 

GRID SIZE TOO LARGE FOR THIS ADAPTION, MINIMUM DIMENSIONS REQUIRED: I D=n,, 
JD=n 2 , KD=/j, (Critical) 

Increase array dimensions to size suggested. Note that these values are appropriate for this set of 
input-control parameters only and may need to be changed for other adaptions of the same grid. For 
multiple grids, this message may be preceded by additional information. 

IMPORT CARD EXPECTED, NOT FOUND (Critical) 

The previous card in the input stream contained export=.t. This must be followed by a set containing 
import =. I. 

INCONSISTENT PLANE AND STEP (Critical) 

A stepping direction has been requested that is not available for the requested plane. For example, 
IKPLANE and JSTEP. 

INPUT FILE SIZES DO NOT MATCH (Critical) 

The grid dimensions on the header records of the grid file and solution file do not match. 

IS Tins FINITE VOLUME? IF SO, SET FV=.TRUE. (Critical) 

The grid file dimensions are one greater than the solution file. Should this be finite volume? 

MAX I TOO LARGE, CHANGE ID TO n, 

MAX J TOO LARGE, CHANGE JD TO n 2 
MAX K TOO LARGE, CHANGE KD TO n, 

CHANGE IMX TO n 4 (All critical) 

These messages occur only for multiple grids, ft implies that one or more of the grids in the file is too 
large for the dimension statement. It need not be the adaption grid. It is possible that even if this is 
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corrected, a subsequent run could indicate a second dimension change to accommodate data swapping for 
the adaption grid. 

NDIM TOO SMALL FOR 2-D TO 3-D TRANSFORMATION (Critical) 

Increase NDIM dimension (e.g., from 4 to 5) so internal transformation can be made. 

NO CONVERGENCE ALONG INITIAL LINE, ERRMIN=fl / (Warning) 

The initial line is a 1-D adaption only. This is rarely a catastrophic error, especially if a l is small; 
however, the adaption may not be completely satisfactory. The only control parameters that affect the 
initial line are RDSMAX, RDSMTN, and NEDGE. 

NO CONVERGENCE ON LINE j AND PLANE k , ERR=a 2 (Warning) 

This message is only a warning and adaption continues. Even many of these messages may of be no 
concern as long as a 2 is small. If adaption is successfully completed, check the new mesh to see if it is 
acceptable. 

NO OUTPUT FILES (Warning) 

For whatever reason, no files have been output. (Is save~.f on all passes? Does Snamel start in column 

2?) 

NO POINTS ADDED (Warning) 

Inconsistency in parameters requesting adding points. 

NO VALUE FOUND FOR QVAL ALONG INITIAL LINE (Critical) 

Outer boundary movement has been requested along a constant contour line. This value has not been 
found on the first grid line. 

NUMBER OF POINTS INCREASED FROM n f TO n 2 (Informational) 

If the ADD option has been input, this message gives the new grid dimensions. 

NUMBER OF POINTS DECREASED FROM n, TO n 2 (Informational) 

If the SUB option has been input, this message gives the new grid dimensions. 

OUTPUT FILES ON UNITS n { AND n 2 (Informational) 

The final grid file will be found on unit n x and the final Q file is written onto unit a? 2 . Note that 
intermediary flies may also have been created. For the simplest single-grid, single-adaption case, default is 
fort. 10 and fort. 1 1. 

PLANE /*,, GRID n 2 COPIED TO PLANE n 3 , GRID n 4 (Informational) 

Indicates the successful transfer of a surface from one grid to another in a multiple grid file. 

POINT(S) REMOVED, NUMBER OF POINTS NOW n (Informational) 

REMOVE option has been invoked. 

^ IS NON MONOTONIC ON LINE j AND PLANE A (Critical) 

This message will terminate the program. It indicates that the values of s, at the completion of the 
iteration on line j are not monotonically increasing, thus implying crossover of points. Since this is 
unacceptable, the program outputs the data. It is then possible to view the plots and re-evaluate the control 
parameters. 

SOLUTION FILE IS F-V: CHECK YOUR BOUNDARY CELLS! (Informational) 

Data stored into the ghost or boundary cells may be incorrect. 

SUB OPTION PRODUCES TOO FEW POINTS FOR ADAPTION (Informational) 

Fewer than 10 points remain in the adaption direction. 

TOO FEW POINTS FOR ADAPTION WITH NEDGE=1 (Critical) 

There are fewer than 10 points along the adaption line. This is not appropriate, especially if NEDGE 
is set. 

TOO MANY MULTIGRIDS FOR NM DIMENSION (Critical) 

The parameter statement at the beginning of each subroutine contains NM, the number of zones in 
a multigrid file. The default value is 10. 
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WARNING: DIRECTIONAL SIGNS INDICATE COORDINATE SYSTEM MAY NOT BE RIGHT- 
HANDED (Critical or informational) 

The calculation of body normals assumes a right-handed coordinate system. In most cases, this 
message indicates the grid should be re-ordered. 

WARNING: ZERO CELL AT (i,j,k), ‘SAGE WILL TRY TO CORRECT!’ 

Routines SHOCKLE and BLCLUST detect zero cells during the reclustering process. This message is 
passed on to the user as a warning. Note that the (i,j,k) are in the transformed coordinate system. The code 
will average two surrounding cells. 


2.6 Outline of Each Subroutine 

The MAIN routine is a driver routine whose task is to call the relevant subroutines. A loop (using 
NADS) is set to provide for multiple adaption passes. The first routine to be called is INITIAL, which reads 
and organizes the data. A loop is then set up for the adaption of each plane, and the constant planar 
variables are computed. Finally, a loop is set up for each line in the plane in which all the coefficients of 
the adaption equation are computed, followed by the solution process. When each pass is complete, the 
OUTPUT routine is called. 

Along with the MAIN routine, the SAGE code consists of the following subroutines, listed here (as 
they are in the source code listing) in alphabetical order. Arguments shown in boldface are computed within 
the subroutine. 

ADDPTS(TER) 

This routine is called by INITIAL if ADD^O. Extra points (depending on ADD) are inserted between 
every grid point within the range LSTADD to LENDADD, by linear interpolation. The error flag is set if 
the additional points exceed the defined dimensions. 

ADDV{C OSX 1 ,COSYl,COSZ 1 ,A 1 ,COSX2,COSY2,COSZ2,A2,COSX,COSY,COSZ) 

This is a utility routine. The direction cosines of two unit vectors (COSXl,COSYl,COSZl) and 
(COSX2,COSY2,COSZ2) are input arguments. The routine computes the direction cosines 
(COSX,COSY,COSZ) of the unit vector that represents the sum of the input vectors, proportioned by the 
coefficients A1 and A2. 

BLCLUST 

This routine is called when the user requests a reclustering of points using the Vinokur algorithm but with no 
outer boundary movement (RECLUST=n>J). Frequently this is used to transform an inviscid grid into a 
viscous grid. New values of s t are computed, and the (x,y,z) and Q values are interpolated. If only a subset of 
the domain is reclustered ( n<imax-l ) the mesh spacing at the boundary is defined by the original spacing 
and not the input value of DSE. 

BLOCK{ J,K) 

Initially, all grid coordinates are stored in the input X, Y, and Z arrays. In this routine, a block of data 
around the current j line is stored in arrays XJ, YJ, ZJ to prevent the original grid from being overwritten 
during interim calculations. Only the converged adaption line is replaced into the original grid arrays. The 
XJ, YJ, and ZJ are made up of all / points on the current j line and all / points on the j-1 and j+I lines in 
the k -l, k and k+I planes, giving (n t , 3,3) dimensioned arrays. At boundaries, nonexistent data points are 
filled with 999. Most calculations within the code, but especially the computation of the normal vectors, 
are performed on this block. 

CLUST2(X, DX0,DX I ,JMAX) 

This computes the new s, (X) based on the algorithm defined by Vinokur (1983). Also includes ASINHF 
and ASINF function routines. 

CROSSV(XT, YT,ZT,XT 1 , YT 1 ,ZT 1 ,DST,COSV, A AP,D AP,ICROSS,J) 

Appendix II, section 1.11.2, contains the analysis used to develop this routine. COSV is the array 
containing the direction cosines of the vector v, defined in that appendix. (XT,YT,ZT) are the 
coordinates of a j line, (XTI,YT!,ZT1) are the coordinates of a juxtaposed line, and DST is A.v along j. As 
an example, this routine is used to compute s-s' (i.e., AA') and DA', the distance between 5 ' and the 
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corresponding node at (ij-l t k), as shown in Fig. 4. ICROSS is the array containing / and indicates the 
intersecting segment for each COSV. 

CSPLIN( NT,S,V,SPF) 

The cubic spline coefficients, SPF, are computed for NT points. S is the streamwise location of the 
function V. These coefficients are used by the routine SPEVAL if GEOM=.true. or INTER=4. 

DETERM{ A I,B1,C1 , A2,B2,C2,A3,B3,C3,DET) 

This routine computes the determinant of the three vectors whose direction cosines are given in the 
argument list. 

DLENG( JL,K) 

When NEDGE is set, the tension parameter co is amended at the edges (i.e., at 1ST and IEND) to improve 
edge-boundary spacing. This routine computes DLENGS and DLENGE, which are used in the edge co 
calculation (by the routine WTEDGE ). The values of DLENGS and DLENGE depend on whether the grid is 
defined outside of the adaption domain. JL indicates which j line is needed. 

EDGE MG (VAR) 

This is a utility routine. For various reasons, the values of some variables at the 1ST and/or TEND edges are 
overridden. To blend these different values into the calculations, this routine will merge the new values 
(given at the two boundaries of VAR) into the next three (or MG1, MG2) grid locations of VAR. 

FBAR( J,K) 

This routine is called once for every j line. The As and the gradients dqjds are computed at the input grid 
points and stored in FQ. If duplicate points are found (i.e., 3 a As, = 0) a message is printed and the 
program terminates. If GEOM-.true., the wall gradients dg/ds are also computed by calling WALLS , and 
stored in FG. In addition, the coefficient C is computed and FG is added to FQ and stored in F. The routine 

INTF is called to interpolate the value of F at these new grid points and to compute the normalized form of 
F, i.e., FB. The routine GETB is called to find the value of B for this j line. For lines other than the first, 
initial guesses for the s distribution are extrapolated from the converged s t along the j-1 line. Since these 
do not correspond to the input points, new local values of x, y, and z are interpolated by calling PROPS. 
FILTERS AR,NIPTS,NFTLT) 

This routine smoothes the parameter contained in VAR by adding a second derivative term, 
y, . = .75v, +. 125(v /+ , - v-_, ). NFILT is the number of smoothing passes (default=2). The parameters 

smoothed are — /{dg/ds) and co,-; they are returned to the calling routine via VAR. 

FVORG{ IND) 

If IND-1, this routine interpolates the cell-centered q values onto the grid points to enable adaption to 
take place in the normal manner. Values at the first and second boundary points are set equal. If IND=2, the 
q values are re-interpolated to the new cell centers. 

GETB(J,K) 

This routine computes the value of B used to evaluate co. B is found by an iterative process and is said to 
converge when the minimum requested As equals the computed minimum As. The analysis for this routine 
is given in Appendix I. 

GETWT 

This routine computes co,, the modifier of co that is applied when any computed As lies outside the 
requested range of RDSMAX and RDSMTN. 

HEADIO( ICJER) 

This routine reads (IC=1) and writes (IC=2) the header records for multiple-grid files. It also tests to see if 
the grid dimensions are too large for the programmed dimensions and whether the NM dimension 
parameter is greater than or equal to the number of zones. 

INLTLA Z,(NOMORE,MT CH) 

This routine sets all the default values and reads the input parameter file. If EXPORT is true for a multiple- 
grid file, MTCFI is set to .true, and control is handed to the MATCLI routine. In all other cases, the 
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appropriate grid and function file input routines are called (either READAT or READMULT). When 
necessary, the grid points and corresponding flow-field data are rearranged to correspond to the data 
organization assumed by the analysis. This routine also controls the addition or removal of points and 
initiates any boundary movement. If no more adaptions are requested, NOMORE is set and control is 
returned to the main routine. 

INTF( F 1 ,F2,S 1 ,SMID,NPTS) 

This routine interpolates to find F2 (new F) at the new s, (SI), based on the input values of FI (current F) 
and the midpoints (SMTD) of the input s array. The interpolation routine chosen is based on INTER. 
INTX}'ZQO,KJ 1 ,K 1 ,SS,SN, QJ) 

This routine interpolates for X, Y, Z, and Q at the new SN, given the corresponding values at SS. The new 
(x,y,z) coordinates are stored in the block of data defined by XJ, YJ, and ZJ. The q data are stored in QJ. 
The appropriate interpolation routine is called based on INTER. 

L4GC0F(SNEW,SARR,NPTS,M,P1,P2,P3) 

This routine computes the Lagrange coefficients PI, P2, and P3 for a point SNEW, with respect to the 
input s array, SARR. These are used by the calling routine to interpolate for the variable at SNEW. First or 
second order is available; the choice depends on INTER. M is the associated index for P I and will reflect a 
forward or backward interpolation, depending on the location of SNEW within the interval. 

LINE 1 

This routine solves for the adapted values of s ( . along the initial line j = j st on the initial plane k sl . Since 
both torsion terms are zero on this line, the As t are computed from the 1-D approach. 

MARCHJ( K) 

If the final adapted line within plane k is internal to the physical end-boundary line, this routine will 
redistribute the points on the remaining j lines, based on the distribution along the j enJ line. The routine is 
called for each plane. This action is performed only on request by the user and is not an adaption to the 
flow field. However, it will prevent the discontinuity between the last adapted line and the remaining non- 
adapted grid lines. 

MARC UK 

The function of this routine is similar to MARCH.! but it is called on user request if the final adaption plane 
k enJ is internal to the physical end plane k max . It will redistribute points on all lines on the remaining planes 
to be proportional to the corresponding line on the last adapted plane. 

MA 7T//(NOMOREJER) 

This routine is called as soon as an ‘export’ card is read. If necessary, the output Files are rewound 
(REWND)\ and then the header records are read ( HEADIO ). The export grid is read into core ( MULTIO ) and 
the requested plane stored ( STOREX ). Now, the next set of input parameters is read (which should be an 
‘import’ set). The input files are closed, rewound, and reread up to the import grid (READMULT). The 
stored export plane is copied into the import plane (STORIKI) and all the files are written to the output 
units (WRITMULT). 

MGWALLS(. J) 

This routine is called when GEOM=.true. It computes the coefficient of the geometry function, C g (FGW), 
based on the local aspect ratio. 

MULTIO(m 1 ,IN2,TOUT 1 ,IOUT2) 

This routine reads and immediately writes multiple grid files. INI and IN2 give the range of grids to read in 
(e.g., if MGRID=3, READMULT will initially call MULTIO with INI = 1 and IN2=3). IOUT1 and IOUT2 
give the range of grids to output (in the same example, READMULT will set IOUTI = l and IOUT2=2, since 
grid 3 will be written after adaption). The header records are handled separately in HEADIO. 

NOADAPT(J,K) 

This routine updates variables and/or places them in appropriate arrays when no adaption is to be 
performed on the current j line, but adaption is to be performed on the next line. When stepping to the 
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next line, the code expects certain variables to be available at /'-/. This scenario occurs, for example, in 
merging situations (MGSTEPS>0) and for common lines (LNSING>0). 

NORM(V 1 ,NPTS) 

The function FI is normalized as (FI - FI „„„)/( FI , Har - FI If maximum and minimum values are equal 
or nearly equal, the normalized variable is set to 0(10 s ). 

M?/?MFF(IP,JP,KP,TNDPL,PA,PB,PC,SING) 

This routine finds the vector at the point (TP,JP,KP) normal to the plane defined by INDPL. Although 
three direction planes exist through the point, only two are needed by the calling routines: INDPL=1 is the 
(i,k) plane and INDPL=2 is the (i,j) plane. The analysis to describe this routine is given in Appendix III. 
Four normals are computed and the average is found, with the direction cosines returned in (PA,PB,PC). 
SING is set to 1 if the normal does not exist. 

OUTPUT 

OUTPUT is called at the conclusion of each adaption set (NADS). Remaining planes are proportioned if 
requested by MARCHK, and then the data files are returned to their original order to conform with the input 
mesh structure. Either WRITMULT or WRITOUT is called to output the grid and flow-field files. This routine 
is also called if 5 becomes non-monotonic (OK=. false.). In this case, the new mesh points that have been 
computed are output to help the user choose more appropriate control parameters. 

PROPS(S,K) 

After a new solution of s, has been obtained on a line /', the code stores this data in line/'-/ and steps to the 
next line. This routine proportions the new s ( throughout the non-adapted regions of the block of data 
defined by XJ,YJ,ZJ. Essentially, this provides a first guess for the current/ line and also produces smoother 
planes for the vector normal calculations (see Appendix III). 

PURPLE{ A 1 ,A2,A3,B 1 ,B2,B3,V1,V2,V3,NFL) 

A and B are direction cosines of two vectors defining an enclosed plane. This routine takes their cross- 
product and normalizes the result to give the direction cosines (VI,V2,V3) of the unit normal to the 
enclosed plane. NFL is the direction sign of the normal. 

mS7//F(DSN,NSM,NSM 1 ) 

This routine controls the outer boundary movement for all /' at each k plane. It is called if MVBOUND^ 0. 
It first calls SHOCKLE to find the new s max . If d n . =999. the original wall spacings are saved for each /'. The 
location of (x,y,z) at the new s max are computed, and then the new array of s t is obtained, either as 
proportional to the old s- or by calling CLUST2. Finally, all (x,y,z) and q are interpolated at the new s- r 

READ AT 

This routine reads in the single grid and function files in PLOT3D binary format. Input datasets may be in 
2-D or 3-D form. The size of IMAX, JMAX, and KMAX on the header record is checked (by SIZING) 
before the data is read. If NOQ is true, the Q array is filled with 1.0. 

READMULT{WR) 

This is the read routine for multiple-grid files. Since some writing of files also occurs, the output unit 
numbers are verified. The header records are read and copied to the output files ( HEADIO ). Then MULTIO 
is called to read the requested grid into core and to write all preceding grid data onto the output files. 
Finally, if NOQ is true, the Q array is filled with 1.0. 

REWND 

This routine rewinds the output files and opens them as input files. 

SETUPJ(I,K) 

This routine computes the direction cosines of the vectors u, b, and e used to evaluate the torsion vector 
7 . These vectors are associated with the (i,j) points within the constant k plane. 

SETUPK( J,K) 

This routine is similar to SETUPJ, but the direction cosines of FT, b* , and e* are computed to evaluate the 
plane torsion vector, 7* . These vectors are associated with the (i,k) points within the constant / plane. 
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StfOCKZ£(SNMAX,DSN,NSM,NSM 1 ) 

This routine is called when the outer boundary is moved, and it computes all the new s mox for each (j,k), 
based on shock location or other user request. Smoothing also takes place, based on input NSM and NSM1. 
SINGPLN 

When PLSING is set, SINGPLN stores the result of the first adapted line into every line of the same plane. 
SIZING(IER) 

The parameter statement at the beginning of each subroutine presets the dimensions (ID,JD,KD) of the grid 
and q files. This routine compares the input grid dimensions (TMAX,JMAX,KMAX) or 
{IM(NGRID),JM(NGRID,KM(NGRID)} to these preset values. If insufficient space has been allocated, the 
minimum possible values of ID, JD, and KD are computed for this adaption to proceed. (These values are 
not obvious since space must be allocated for data swapping.) IMX, the maximum of ID, JD, and KD is also 
evaluated. A message is then sent to the user to recompile SAGE with the suggested dimensions. Finally, 
IER is set to 1 to inform the INITIAL routine to terminate the code. 

SOLUT(I, K) 

By the time SOLUT is called, all variables have been computed that are needed to obtain the new 
distribution of s r The coefficients of s, (see Eq. (8) in the first section of this report) are set up in a 
tridiagonal matrix and solved for s,. Interpolated values of w, are found at these new values of s, and 

iterations are performed until is small or too many iterations have been performed. In both 

cases, a check is made to see if all the s t are monotonic. If so, the program continues; if not, the flag OK is 
set to false, causing the program to output the current grid and terminate. 
SP£UIL(NT,S,V,SPF,SI,VI,VPI,VPPI) 

This is a cubic spline interpolation routine used when TNTER=4 or if GEOM is true. It uses the coefficients 
(SPF) computed in CSPLIN. In addition to interpolating for the variable V at SI, the corresponding first 
(VPI) and second (VPPI) derivatives are also evaluated. 

STOREX 

This routine stores the data from the export plane into XP, YP, ZP, and QP. 

STORIM 

The export data stored by STOREX in XP, YP, ZP, and QP are copied to the import plane. 

SUBPTS 

When SUB = n ^0, n points are deleted between each retained mesh point. If LSTSUB and/or LENDSUB 
are nonzero, the range of the deletion is restricted. Deletion occurs only in the adaption direction and does 
not decrease the number of stepping lines or planes. 

SWAP IN V 

This routine is called if k st > k enJ , j sl > j rml , or i sl > i em , . The order of (i,j,k) in the x, y, and z input matrices is 
reorganized to ensure that internal computations have monotonically increasing indices. The flag for the 
handedness of the coordinate system is amended accordingly. 

SIUIP.YT^RSWAP) 

Since the internal computation assumes that j is the stepping direction within the plane and that k is the 
stepping direction of planes, this routine is called to interchange x, y, z, and q when the input requests 
alternative stepping directions. RSWAP is a flag that states whether this data exchange is at the start of the 
computations or is the reverse process required for output. This routine is lengthy due to minimizing 
storage requirements. 

SWAP2D(\Q) 

TWOD=.true. indicates that the input datasets are formatted in two dimensions. This routine reorganizes 
the 2-D plane to appear in the code as a 3-D surface. Every' z is given the value of zero, each Q(4) is moved 
to Q(5), and Q(4) is zeroed. 10 indicates whether the translation is from 2-D to 3-D, which occurs on 
input, or from 3-D back to 2-D for output. 
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r07?C0F(L,JK,JKST,JKEND,MGN0S,MARCHJK) 

TORCOF is called twice in each loop, once with all the arguments representing the k plane passing through 
(ij, k) and once with the arguments representing the j plane also passing through (ij,k). This routine amends 
the coefficients (C,,X) etc. of the torsion vectors i and t* based on the current line location. For example, 
C t is decreased when leaving or approaching a boundary to emphasize orthogonality. 

TORSION( J,K) 

This routine first chooses the appropriate b and adds to u to obtain h. The torsion vector t is then 
obtained by adding h and e . The between-plane torsion vector t* is computed in a similar manner. The 
routine CROSSV is then called to find the intersection of the torsion vectors with the j line from which 
both s' and s* can now be evaluated. Finally, a check is made to ensure that s' and s* monotonically 
increase. If they do not, the code attempts a correction, but any major problems will cause the code to 
terminate in the SOLUT routine. 

UNITV(X 1 ,Y 1 ,Z 1 ,X2,Y2,Z2,DIRCX,DIRC Y,DIRCZ) 

This is a utility routine that finds the unit vector from (X1,Y1,Z1) to (X2,Y2,Z2). DTRCX, DTRCY, and 
DIRCZ are the direction cosines of this vector. 

UPDATE^, K) 

This is the last routine called in the iteration loop for the current j. Newly adapted values of s i have been 

found. The values of x, y, z y and q at this new distribution are interpolated (INTXYZQ) and replaced into the 
matrices containing the physical mesh. 

EM£/?G£(DIRV,LST,LEND) 

This routine performs the same function as EDGEMG , but with a vector quantity (DIRV) in place of a 
scalar value. LST and LEND indicate which value of DIRV must be merged into the next three (or MG1, 
MG2) points. 

WALLS{ JW,K) 

Along line j (JW), the geometry gradient, as defined by the radius of curvature, is computed for each grid 
segment. Normally, j w ill equal j st or j end . However, for cases when geometry is the only adaption variable 
(QFUN=. false.), WALLS is called for every j line. 

WR1TMULT 

This is the write routine for multiple-grid files. The current adapted grid is now output and any subsequent 
grids are read and written by calling MULTIO . If the grid size has changed, the header record is updated; the 
output files are closed and opened as input files; the header record is corrected and all data are copied to the 
new output files. 

WRITOUT 

This routine writes the single-grid 2-D or 3-D adapted grid and interpolated function file on units NITG and 
NITQ. If the finite-volume option is set, IMAX, JMAX, and KMAX are one less on the q file. 

WTEDGE( J,K) 

This routine is called by MAIN when NEDGE modification is requested. The edge values of the tension 
parameter at the next j line are a function of the average coA s along the just-completed adapted line. This 
routine calls DLENG to obtain the appropriate value of edge As on the next line, computes the average 
COA s on this line, and evaluates WDS and WDE to be used in routines LINE1 and SOLUT. 

2.7 Nomenclature 

The following is a list of variables used in the SAGE analysis. When applicable, the corresponding 
FORTRAN name used in the code is shown in boldface. 

A, B constants used to compute co, (A,B) 

A v aspect ratio, used to control C7,,FGASP 

C g coefficient of f(g) in co calculation, FGW 
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coefficient of f(q) in 03 calculation, FQW 
input proportion coefficient for torsion, CT(1) 
input proportion for torsion between planes, CT(2) 

orthogonal vector to j line within k plane; direction cosines of b , COSB 

orthogonal vector to j line within / plane, COSBK 
modified values of C n C* , CTM(l), CTM(2) 

straightness vector, COSD 

average straightness vector within k plane, COSE 

average straightness vector within j plane, from k — 2 — > k — 1 , COSEK 
gradient of q (and g if necessary), F 

minimum and maximum / used to normalize f FMIN, FMAX 
normalized function of /, FB 
geometry function 

subscript indicating the current node in adaption direction, I 
total number of points in i,j and k directions of input grid file, 

IMAX, JMAX, KMAX 

start of adaption domain in i,j and k directions, 1ST, JST, KST 

indices indicating end of adaption domain, TEND, JEND, KEND 

subscript of the current stepping line, J 

subscript of the current adaption plane, K 

torsion-related constant 

local subscript relating to node /, L 

number of stepping lines before full adaption, MGSTEPS 

the plane equivalent of MGPLS 

flag for edge control, NEDGE 

number of points in the adaption line, NTPTS 

merging coefficient for lines f(m g ) 9 CNM(l) 

merging coefficient for planes /(m*) 3 CNM(2) 
orthogonality vector within plane, COSN 

orthogonality vector between planes, also stored in COSN 

input flow-field variable (p t pu,pv,pw,e), Q 

radius of curvature for geometry 7 function, FGS 

streamwise length, SS or SN 

vector representing s 9 SSX, SSY, SSZ 

maximum value of s on line j, SMAX 

value of streamwise length used for torsion within planes, SP 
value of streamwise length used for torsion between planes, SPP 
DS 

requested minimum and maximum grid spacing, DSM1N, DSMAX 
computed minimum and maximum grid spacing 
torsion force 

within-plane torsion vector, COST 
torsion vector between planes, also COST 
proportion of b and Ft used to compute n,TN(l) 
proportion of b* and u* used to compute Ti\ TN(2) 
vector normal to j-1 line in k plane, COSU 


36 



II 

x,y,z 

X 

X* 

0 

x 

(0 

CO, 

V 


vector normal to j line on k-1 plane, COSUK 

input grid mesh, (X,Y,Z) globally and (XJ,YJ,ZJ) locally 

input magnitude of torsion control parameter, CLAM(l) 

input magnitude of plane torsion control parameter, CLAM(2) 

angle for torsion computation 

witliin-plane torsion-related parameter, TAU 

tension spring force, WEIGHT 

computed weighting on co, WT 

between-plane torsion-related parameter, TAUPL 


2.8 List of Major Variables 


This section contains the list of variables (in alphabetical order) used in the SAGE code. Local variables 
that contain only intermediary values are not listed. The format is: 

Variable name(dimension) /r, / r 2 /brief description 

where describes what type of variable — input, local, parameter, or common block name 

and r 2 lists routine(s) where the variable is initialized. 


A 

AA(IMX) 

AAPfIMX) 

ACT 

ADD 

ALPHA 

AMACH(TMX) 

B 

BB(IMX) 

BCONV 

BJ1 

BKI 


/COM2/INITIAL/ A used to compute 0) 

/local/SOLUT/ coefficient of s,_| in solution matrix 
/argu nient/CROSSV/ s — s' or s — s* 

/local/TORSION/ final modified C, 

/COM 1 9/input/ if set, add grid points 
/COM1 2/input/ information only for PLOT3D 
/local/FBAR/ computed Mach number 

/COM2/GETB/ B used to compute (0 
/local/SOLUT/ coefficient of s t in solution matrix 
/local/GETB/ convergence criteria for B iteration 
/local/GETB/ value of B along/-/ line in the k plane 
/COM15/GETB/ value of B along/ line in the k-1 plane 


CC(TMX) 

CLAM(2) 

CLAMW(2) 

CNM(2) 

CONV 

COSB(lMX,3) 

COSBK(TMX,3) 

COSD(IMX,3) 

COSE(IMX,3) 

COSEK(IMX,3) 

COST(TMX,3) 

COSU(IMX,3) 

COSUK(TMX,3) 

CT(2) 

CTM(2) 


/local/SOLUT/ coefficient of s, +l in solution matrix 
/COM 10/input/ X and X , magnitude of torsion terms 
/COM 1 0/TORCOF/ modified CLAM 

/COM I O/TORCOF/ X modifiers for MGSTEPS* 0 and MGPLS * 0 
/COM 1 5/INITIAL/ general convergence criteria 
/COM7/TORSION/ direction cosines of b , in (i.j) plane 
/COM20/TORSION/ direction cosines of b * , in (i,k) plane 
/local/SETUPJ,SETUPK/ temporary straightness vector, d 
/COM7/SETUPJ/ straightness vector, e , in (ij) plane 
/COM20/SETUPK/ straightness vector, e* , in (i,kj plane 
/local/TORSION/ torsion vector, t or 1* 

/COM7/SETUPJ/ u , normal to j-I line, in k plane 
/COM20/SETUPK/ «*, normal to / line in k-1 plane 
/COM 1 0/input/ C, and C* , directions for torsion vectors 
/COM I O/TORCOF/ modified C,and C] 


DAP(IMX) 

DAPPL(IMX) 

DET 

DLENGE 


/COM9/CROSSV/ DA' for 5 ' calculation 
/COM9/CROSSV/ DA' for U calculation 
/argument/DETERM/ value of three-order determinant 
/COM6/DLENG/ As computed for end-edge modification 
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DLENGS 

DMINSDB 

DSN 

DS(TMX) 

DSE 

DSMAX 

DSMIN 

DSW 

DSW1,DSW2 

EXPORT 

F(IMX) 

FB(IMX) 

FF(IMX) 

FG(IMX) 

FGASP 

FGS(IMX) 

FGW 

FORMI 

FORMO 

FQ(TMX) 

FQW 

FSMACH 

FV 

GEOM 

ICROSS 

ID 

IEND 

IINVERSE 

IJPLANE 

IKPLANE 

IM(NM) 

IMAX 

IMAXQ 

IMPORT 

IMQ(NM) 

IMX 

INDQ 

INTER 

IQ(NDIM+3) 

IS,IE 

1ST 

I STEP 

JD 

JEND 

JINVERSE 

JKPLANE 

JM(NM) 

JMAX 

JMAXQ 

JMQ(NM) 

JS,JE 


/COM6/DLENG/ As computed for start-edge modification 
/local/GETB/ dmin(As)/dB 

/argument/input,PUSHIT/ parameter for outer boundary location 

/C0M3/FBAR,LINE1 ,SOLUT/ As 

/COM 1 8/input/ outer edge spacing (when RECLUST * 0) 

/COM6/FBAR/ As max from RDSMAX 

/COM6/FBAR/ A s min from RDSMIN 

/COM 1 8/input/ wall spacing for Vinokur algorithm 

/COM 1 8/input/ first and last wall spacing when DSW=999 

/COM21 /input/ if true, store data for transferring between multigrids 

/COM2/FBAR/ flow gradient, / = dq/ds , at input s 
/COM2/INTF/ / , normalized / at current s 
/local/SOLUT/ coefficient of RHS of solution matrix 
/COM 1 7/FBAR/ normalized dg/ds 
/local/MGWALLS/ function of aspect ratio 
/COM 17/WALLS/ dg/ds at SMSS 

/COM1 7/MGWALLS/ coefficient of FG for to calculation 
/COM24/input/ true if ASCII (formatted) input files 
/COM24/input/ true if ASCII (formatted) output files 
/COM 1 7/FBAR/ normalized dq/ds 

/COM1 7/MGWALLS/ coefficient of FQ for to calculation 
/COM 1 2/input/ information only, for PLOT3D 
/COM 1 5/input/ flag to indicate finite-volume q file 

/COM1 1 /input/ request for geometry fiinction 

/argument/CROSSV/ index for interval location of s' or s * 
/dimension/parameter statement/ grid dimension in i direction 
/COM4/input/ last node along i in adaption domain 
/COM 13/INITIAL/ adaption requested in backward / steps 
/COM 14/input/ true if ( ij ) is adaption plane 
/COM 14/input/ true if (i,k) is adaption plane 
/COM2 1 /input/ IMAX for each grid in multigrid file 
/COM4/input/ number of points in physical domain, / direction 
/COM24/input/ IMAX for q file 

/COM2 1/input/ if true, describes location to place EXPORT data 
/COM24/input/ IMAX for each grid in the multigrid q file 
/dimension/parameter statement/ maximum value of ID,JD,KD 
/COM5/input/ index for q for adaption variable 
/COM1 1/input/ order of interpolation, 2, 3, or 4 
/COM5/input/ related to INDQ, combines q adaption function 
/COM23/input/ i range of export/import transfer plane 
/COM4/input/ first node along / in adaption domain 
/COM 13/input/ true for stepping in / direction within the plane 

/dimension/parameter statement/ grid dimension iny direction 
/COM4/input/ last node in j adaption domain 
/COM 13/INITIAL/ adaption requested in backward j steps 
/COM 14/input/ true if (j.k ) is adaption plane 
/COM21 /input/ JMAX for each grid in multigrid file 
/COM4/input/ number of points in physical domain,./ direction 
/COM24/input/ JMAX for q file 

/COM24/input/ JMAX for each grid in the multigrid q file 
/COM23/input/ j range of export/import transfer plane 
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JST 

JSTEP 

KD 

KEND 

KINVERSE 

KM(NM) 

KMAX 

KMAXQ 

KMQ(NM) 

KS,KE 

KST 

KSTEP 


/C0M4/input/ first node in j adaption domain 
/COM13/input/ true for stepping iny direction within the plane 

/dimension/parameter statement/ grid dimension in k direction 
/COM4/input/ last node in k adaption domain 
/COM 13/INITIAL/ adaption requested in backward A: steps 
/COM21 /input/ KMAX for each grid in multigrid file 
/COM4/input/ number of points in physical domain, k direction 
/COM24/input/ KMAX for q file 

/COM24/input/ KMAX for each grid in the multigrid q file 
/COM23/input/ k range of export/import transfer plane 
/COM4/input/ first node in A: adaption domain 
/COM 13/ true for stepping in A: direction within the plane 


LEND ADD 

LENDSUB 

LNSING 

LSTADD 

LSTSUB 


/COM 19/input/ start of add points range 

/COM 19/input/ start of delete points range 

/COM! 5/input/ flag to indicate a common line for each plane 

/COM 19/input/ end of add points range 

/COM 19/input/ end of delete points range 


MARCH 

MARCHPL 

MAXTTS 

MG1 

MG2 

MGCT(NM) 

MGPLS 

MGRID 

MGSTEPS 

MPLANE 

MSI,MSJ,MSK 

MTCH 

MVBOUND 

MXFG 


/COM 14/input/ true to interpolate up to physical boundary 
/COM 14/input/ true to interpolate to final plane 
/COM 15/INITIAL/ maximum number of iterations for convergence 
/COM9/INITIAL/ number of merging points for NEDGE at i xl 
/COM9/INITIAL/ number of merging points for NEDGE at i end 

/COM 16/INITIAL/ counter on number of adaptions per grid 
/COM1 1/input/ number of planes before full adaption within plane 
/COM12/input/ multiple grid number for adaption 
/COM1 I/input/ number of lines before full adaption within plane 
/COM23/input/ plane number for export or import plane 
/COM21/SIZING/ contains overall maximum size of multiple grids 
/argument/INITIAL/ Hag to indicate export and import has occurred 
/COM1 8/input/ invokes outer boundary movement 
/COM 17/MG WALLS/ location of maximum dg/ds 


NA 

NADS 

NDIM 

NEDGE 

NEILT 

NFLAG 

NGRTD 

NIPTS 

NITGI 

NITGO 

NITQI 

NITQO 

NM 

NOMORE 

NOQ 

NOSESM 

NOUP 

NSM/NSMI 


/dimension/parameter statement/ maximum number of adaption passes 
/COM4/MAIN/ index on number of adaptions 

/dimension/parameter statement/ no. of q variables, normally set to 5 
/COM9/input/ initiates side-edge boundary override 
/COM1 1 /input/ number of passes for smoothing data 
/COM 15/INITIAL/ indicates direction sign for normal vectors 
/COM2 1/input/ number of grids in multigrid file 
/COM4/INITIAL/ total number of / points in computation domain 
/COM 1 2/INITIAL, WRITOUT.../ unit number for input grid file 
/COM12/INITIAL,WRITOUT.../ unit number for output grid file 
/COM 12/INITIAL, WRITOUT.../ unit number for input q file 
/COM 12/INITIAL, WRITOUT.../ unit number for output q fde 
/dimension/parameter statement/ maximum number of multigrids 
/argument/INITIAL/ no more input datasets remain; end program 
/COM1 1 /input/ false to indicate no q file will be input 
/COM 15/input/ points to smooth nose kink after moving boundary 
/COM1 I/input/ prevents adaption, i.e., no update 
/argument/input, PUSHIT/ outer boundary smoothing 


OK 

ORTHE(2) 

ORTHS(2) 


/COM14/SOLUT/ flag to indicate normal termination 

/COM 16/input/ false removes orthogonality at outer boundaries 

/COM 16/input/ false removes orthogonality at initial boundaries 
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P1,P2,P3 

PLSING 

PRES(IMX) 

Q(ID,JD,KD,NDIM) 

QFUN 

QP(IMX,IMX,NDIM) 

RDSMAX 

RDSM1N 

RE 

RECLUST 

REMOVE 

RSWAP 

SAVE 

SING 

SMS(IMX) 

SMSS(IMX) 

SN(IMX) 

SNM(TMX) 

SNMK(IMX) 

SP(IMX) 

SPF(IMX) 

SPPL(TMX) 

SS(IMX) 

SUB 

TAU 

TAUPL 

TIME 

TN(2) 

TNM(2) 

TRAT(IMX) 

TWOD 

WDE 

WDES 

WDS 

WDSS 

WEIGHT(IMX) 

WT(IMX) 

WTSUM 

X(ID,JD,KD) 

XJ(IMX,3,3) 

XP(IMX,IMX) 

Y(ID,JD,KD) 

YJ(IMX,3,3) 

YP(IMX,TMX) 

Z(ID,JD,KD) 

ZJ(IMX,3,3) 

ZP(IMX,IMX) 


/argument/LAGCOF/ coefficients of Lagrange polynomials 
/C0M1 5/input/ flag to define a singular plane 
/local/FBAR/ computed pressure 

/COM 1 /input, UPDATE/ flow- field variables 

/COM1 1 /input/ false to allow geometry function only 

/COM23/STOREX/ q data of stored plane for grid transfer 

/COM5/input/ relative value of As max 

/COM5/input/ relative value of A s min 

/COM 1 2/input/ not used, PLOT3D variable 

/COM 1 8/input/ invokes Vinokur reclustering algorithm 

/COM 18/input/ number of points to remove at outer boundary 

/argument/INITIAL,OUTPUT/ flag for swapping data 

/COM 1 4/input/ flag to suppress output 
/argument/NORMPT/ flag to indicate nonexistent normal 
/COM3/FBAR/ midpoints of As 
/COM I 7/WALLS/ s, for geometry calculation 

/COM3/FBAR/ current s array alongy line 
/COM3/UPDATE/ converged s array alongy'-/ line 
/COM3/FBAR/ converged s array alongy line on k-1 plane 
/COM9/TORSION/ s' , intersection of torsion vector and y line 
/argument/C SPLTN/ cubic spline coefficients 

/COM9/TORSION/ j* , intersection of plane torsion vector and j line 
/COM3/FBAR/ initial value of s alongy line 
/COM 1 9/input/ request to delete points 

/local/SOLUT/ t, torsion parameter within plane 

/local/SOLUT/ \|/, torsion parameter between planes 

/COM 1 2/input/ not used, PLO 131) variable 

/COM 1 0/INITIAL/ proportion of b to « for torsion vectors 

/COM 1 0/TORCQF/ TN modified for boundaries 

/local/FBAR/ computed temperature ratio 

/COM 13/input/ indicates 2-D input files 

/COM2/WTEDGE/ weighting factor used when NEDGE set 
/COM6/WTEDGE/ retained WDE 
/COM2/WTEDGE/ as WDE, but at initial boundary 
/COM6/WTEDGE/ retained WDS 
/COM2/LINEl,SOLUT/ co, tension parameter 
/COM6/GETWT/ co,, correction to co 
/local/GETB,LINE 1 ,WTEDGE/ ^ I/co 

/COM ] /input/ input grid points, / direction 
/COM8/BLOCK,INTXYZQ/ X data stored in block 
/COM22/STOREX/ stored X of export transfer data 

/COM I /input/ input grid points,y direction 
/COM8/BLOCKJNTXYZQ/ Y data stored in block 
/COM22/STOREX/ stored Y of export transfer data 

/COM 1 /input/ input grid points, k direction 
/COM8/BLOCK,INTXYZQ/ Z data stored in block 
/COM22/STOREX/ stored Z of export transfer data 
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3. EXAMPLES 


This section contains many examples to familiarize the user with the adaptive-grid process. Each 
example includes plots of the initial grid and flow-field contours, the input-control-parameter file used to 
adapt the grid, the resultant adapted grid, and a discussion on the choice of control parameters. In addition, 
some cases show the improved flow solution obtained from the flow solver using the adapted grid as input. 
The 2-D examples show the effect of within-plane parameters that are also fundamental to 3-D problems, 
so 3-D users should also study these examples. 

3.1 Two-Dimensional Examples 

The first set of examples is for 2-D problems. Although the code adapts the data as if it were a 
single plane in a 3-D environment, this is imperceptible to the user. When TWOD= true. in the input 
parameter file, SAGE will expect the data to be in the 2-D format of PLOT3D. The index of (1) on the 
CLAM and CT parameters refers to within-plane variables. Index (2) is not needed for these 2-D cases. If 
TWOD=.true. on the first pass, it will be assumed for subsequent passes in the same run. 

Case 1. Flow in a Supersonic Inlet 

Figure 13(a) shows the 101 x79 initial grid (assigned to unit 7) for an inlet flow-field problem. Flow 
is from left to right (defined as the / direction) and a shock emanates from the upper-wall comer, reflecting 
off the lower wall. In addition, an expansion fan originates from the downstream upper-wall comer and 
interacts with the reflected shock. An interim computed solution (input on unit 8), generated by the flow 
solver using this initial grid, is shown as density contours in Fig. 13(b). 



Figure 13. Flow in a supersonic inlet, (a) Initial grid; (b) computed density contours. 

The SAGE code is now run to create a new grid that is more adapted to this solution, i.e., the grid 
points are redistributed with respect to a chosen flow-field variable and output to the file assigned to 
unit 10. SAGE also interpolates the complete flow-field solution onto these new grid points and outputs this 
file to unit 1 1. The updated files will subsequently be input into the flow solver to produce more accurate 
solutions. Several examples of grid adaption are given for this inlet problem to demonstrate the effect of 
varying the control parameters. Each example uses the grid and solution files shown in Fig. 13 as input. 

Example 1: Single pass, stepping in j direction. Until the user is familiar with the basic parameters, 
the first attempt should be an input-control file with no parameters; i.e., all default values are used (with the 
exception of TWOD=.true.): 

Snamel twod=.t. $ 

On completion of the execution, the message to unit 6 states "ADAPTION I COMPLETE". The adapted 
grid using the default input parameters is shown in Fig. 14(a). The grid has been more evenly spaced and 
there is a slight clustering of points around the shock on the lower-wall boundary but the remaining grid 
lines are not adequately adapted to the flow features. This implies that the torsion term (whose magnitude 
is controlled by CLAM(l)) is too large and is overriding the local tension term, preventing the tension 
forces from pulling the points to the high-gradient regions. In addition, the minimum grid spacing is too 
large. Based on this information, a new input-control file was chosen: 
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Snamel twod= . 1. , rdsmin =. 25, clam(l) =.001, nedge = 1 $ 

Only three control parameters are input: the minimum allowable grid spacing, the magnitude of the torsion 
term (an order of magnitude less than the default value), and a request for edge control (which is frequently 
used). The remaining parameters retain their default values, signifying that the full grid will be adapted, 
density is the adaption variable, and stepping is in the j direction. 




Figure 14. Adapted grids, (a) All default parameters; (b) adaption from bottom to top: 

rdsmin=.25, clam(I)=.001, nedge 1 . 

The result of this adaption is shown in Fig. 14(b). It can be seen that, since adaption began along the lower 
surface (Jst=J is the default), the redistributed points cluster around the point of reflection on this line. 
However, points do not become sufficiently clustered at the comer shock and at the start of the expansion 
wave region on the upper surface. This is to be expected, since the adaption on the initial line is a l-D 
solution and will pick up features quite clearly. However, as stepping continues, the features on subsequent 
lines get "dampened" by the torsion (i.e., smoothness) control. 

Example 2: Adaption with backward y'steps. This example maintains the same parameters as Ex. 1, 
except that the upper surface is chosen as the initial adaption line, i.e., 

$ name l twod , rdsm in =.25, clam (I) .001, nedge = IJst = 7 9,jend= 1 S 

Figure 14(c) shows the result of this adaption. In this case, the upper surface is more clearly adapted, and 
the point of reflection on the lower surface is more spread out. Comparing these two examples indicates 
quite clearly that the starting line has a strong effect on the resultant adapted grid, and for some 
applications it should be chosen carefully. • : r 



Figure 14 concluded, (c) Marching from top to bottom; (d) marching in i, from right to left. 

Example 3: Adaption in / direction. This example shows the very different grid generated when the 
adaption is performed by stepping in the i direction ( istep=.true .). Based on the experience obtained from 
the first two examples, adapting from right to left will produce a better grid since there are no gradients on 
line ist=l to adapt to. Hence, the input control file was chosen to contain: 

Snamel twod=. t . , istep=.f, ist=101, iend= l.rdsmin =.25,clam(l) =. 001, nedge=l S 
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Figure 14(d) shows the result of this adaption. Not surprisingly, the reflected shock region on the lower 
surface is not ’’captured” by this adaption, and thus this grid is less suitable than those shown above. This 
demonstrates that the choice of stepping direction is important in producing a desirable grid. 

Example 4: Effect of changing control parameters. This example is actually a collection of examples 
that show the effect of varying the control parameters clam(l), ct(l), rdsmax, and rdsmin. Figure 15(a) 
shows an adapted grid using "baseline” values. These are the same as the default values in the code with the 
exceptions of clam(l)=.001 and adaption proceeds from top to bottom (Le.,/s/=79, jend=l ). 

Snamel twod=. true . , clam(l)=. 001 Jst = 79,jend= 1, nedge - 1 $ 



Figure 15. Effect of control parameters, (a) Baseline: clam(l)=.001 t ct(l)=.5,rdsmax=2.0,rdsmin=.5 ; 

(b) clam(l)=.0I ; (c) clarn(l)=.0001; (d) ct(l)-L0 ; (e) ct(l)=.0 ; (f) rdsmax -4.0, rdsmin-. 25. 

Figures 15(b) through 15(f) are grids adapted by changing just one of the baseline parameters. Figure 15(b) 
shows the result with clam(l)=.0f and it is immediately obvious why this code default value was not 
chosen as the baseline value for these comparison cases: clam(l) is too large to allow any of the flow 
features to be "captured” by the adaption. Figure 15(c) shows the case for clam(l)^.0001 , and this smaller 
value produces a grid with more points clustered around the shocks. Figures 15(d) and 15(e) show the effect 
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of the ct(I) parameter: ct(I)=1.0 in Fig. 15(d), emphasizing straightness, and ct=.0 in Fig. 15(e), 
emphasizing orthogonality. The latter shows how the orthogonality term dampens the adaption for this 
case: we are requesting orthogonality to the already parallel / lines. Finally, the effect of changing the 
minimum and maximum allowable grid spacing is shown in Fig. 15(f), where rdsmax=4. 0 and rdsmin =.2 5 
are used. This mesh size control is only a factor of two different from that in the baseline example (Fig. 
15(a)), but significantly densifies the spacing in the shock regions. 

Example 5: Two-directional adaption. Two-directional adaption is created by adapting in one direction 
(e.g., stepping in j) and then adapting in the other direction (i.e., /), using the first adapted grid as input. 
The resultant grid will depend on the order of the stepping direction. Two passes (i.e., two sets of input 
control parameters) were made to produce the adapted grid shown in Fig. 16(a). If both these passes are 
made during the same execution run (not to be recommended for a first attempt), then additional output 
files will be created on units 12 and 13. Data will be output for both passes, unless the SAVE parameter is 
set to false on the first set of control parameters. The control file for this example contains: 

Snamel twod=.t.,clam(l) =. 001, rdsmax=4. 0, rdsmin 25Jst=79,jend=l, nedge = 1 $ 

S name I istep =. /. , ist =101, iend= 1, clam(l ) =. 001, rdsmin =.25, nedge =1 $ 

The first pass steps in the j direction and uses the same input control parameters as the adapted grid 
shown in Fig. 15(f). The second pass steps in /, and uses the same parameters as Ex. 3. Note that any 
parameters that remain unchanged for the two passes still need to be redefined in the control file, since the 
code restores all default parameters at the conclusion of each pass (with the exception of twod=.true.). The 
control file that produces Fig. 16(b) contains the same parameters, but the adaption order is reversed. That 
is, the first pass steps in the / direction and the second pass in the j direction. The difference between Figs. 
16(a) and 16(b) is obvious and shows that the initial grid completely effects the adaption. Although this 
example shows that both adaptive sequences do adapt the grid quite reasonably, it should be noted that it is 
also possible that one order of adaption will produce a completely unacceptable result, while the reverse is 
quite suitable. 



Figure 16. Two-directional adaption, (a) Marching in j followed by marching in i; 

(b) marching in i followed by marching in j. 

Example 6: Flow solution using adapted grid. Figure 17 demonstrates the improved flow-field 
solution obtained when the flow solver is rerun using an adapted grid as input. The adapted grid shown in 
Fig. 17(a) was obtained by choosing the ‘best’ parameters, based on the experience already described: 

Snamel twod=.true. ,jst = 79,jend= l,rdsmax=4. 0, rdsmin =.25, clam( 1) =. 0001, nedge =1 $ 

This grid and the interpolated flow-field variables were input to the flow solver, producing the flow solution 
(shown as density contours) seen in Fig. 17(b). The improvement in the resolution of the incident and 
reflected shock is obvious when compared to the initial solution shown in Fig. 13(b). 
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Figure 17. Solution using adapted grid as input, (a) ‘ Best ' adaption , rdsmax=4.0 t rdsmin~.25 t 
clam(])=.0001; (b) solution density contours , using ‘ best ’ adapted grid. 

Case 2. Hypersonic BIunt-Body Flow 

The second case contains two examples that demonstrate the effect of the CLAM and CT 
parameters. Figures 18(a) and 18(b) show an initial grid (32x32) and corresponding density flow-field 
contours for a hypersonic blunt body. The j direction marches out from the surface to the free-stream 
boundary, and the i direction marches along the body starting at the lowest point. This is a simple one- 
directional problem with the shock shape aligned with the grid. The outer-side boundary is curved (when 
marching in the / direction) and the default values of clam=.01 and ct=.5 will prevent the adapted grid lines 
from "turning” sufficiently at this edge, giving either a false densing of points at the outer side boundary or, 
possibly, a critical error message. In this situation, two remedies are available. CLAM can be decreased, 
hence de-emphasizing the effect of torsion and thereby allowing the tension force to ’'pull" the nodes 
toward the shock wave. Alternatively, the orthogonality restraint can be removed by setting ct = 1.0. Figure 
18(c) shows the result of setting clam=.001 and retaining all other parameters at their default value. 



Figure 18. Hypersonic blunt body, (a) Initial grid; (b) initial density contours; 

(c) adapted grid. clam-. 001. 
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The adapted grid obtained by setting ct=1.0 is very similar to that shown in Fig. 15(c). The following are 
the two input-control files: 

Snamel twod=. t. , istep=. t. , clam(l )=.001S and Snamel twod=. t. , islep =, t. , cl (I) =1.0 $ 

Both adapted grids show the required result, i.e., the clustering of grid points across the shock. 


Case 3. Blunt-Body Shock Impingement Problem 


The third case introduces the use of the IQ and ORTHE parameters. Figure 19(a) shows the input 
grid of a cowl/lip shock interaction problem. The j direction marches from the body surface to the outer 
free stream and the / direction marches from the lower point of the sphere. The flow-field contours of both 
density (Fig. 19(b)) and Mach number (Fig. 19(c)) are given, since the adaption flow-field parameter is a 
combination of the two. The blunt-body shocks, impinging shock, and shear layers represent a complex 
flow that requires adaption in both directions to adequately capture all the flow features. Figure 19(d) shows 
the adapted grid (output on unit 13) produced by the following two-pass control file: 

Snamel twod=. I. . islep =. I. , isl =91, iend=l, indq=0, iq(l) =2, iq(7) = 1, rdsmin =.25, 
cl=. 7,clant(l)=.005,nedge=2 S 
Snamel indq =0. iq(l) =2, iq(7j— 1 , rdsmin =. 3, or I he ~.f S 

The First pass steps in the i direction, starting from the topmost boundary. The parameter indq=0 
indicates that IQ will be input, and in this case two-thirds of the density function and one-third of the Mach 
function will be combined to become the adaption variable. As explained in case 2, the curvature of the 
outer boundary requires increasing CT. Setting nedge=2 requests that only the side-edge boundary at j=l be 
treated. The second pass introduces the use of the ORTHE variable: stepping occurs from the sphere wall to 
the outer curved boundary, where the default would turn the grid to be normal to this outer boundary. 
ORTHE overrides this and allows the adaption to occur naturally. 



Figure 19. Blunt-body shock impingement problem, (a) Initial grid; (b) initial density contours; 

(c) initial Mach contours; (d) adapted grid. 


Case 4. Hypersonic Inlet (Zonal Adaption) 

Figures 20(a) and 20(b) show the initial grid and density contours for a hypersonic cowl/lip inlet 
problem. This is a more complex case and requires dividing the adaption domain into two zones: the blunt 
body region and the rectangular inlet region. The upper wall is the j=l line and the outgoing channel region 
(on the right side of the diagram) is the i=l line. 
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Figure 20 . Hypersonic inlet , zonal adaption, (a) Initial grid; (b) initial density contours; (c)adapted grid. 


Five adaption passes are required to create the final adapted grid shown in Fig. 20(c). The input 
control file consists of 

Snamel twod=.t. t rdsmin=l. 0,rdsmax=l. 0,save-.f $ 

5 nam e 1 jstep ~.f., rdsmin = 2, nedge = 1, save =./ 5 
Snamel ist=70,save=.f S 

Snamel jst=32Jend=l f lend = 77, rdsmin 25, clam(l)--. 02, nedge - 1, save =/ S 
Snamel iend = 85,jst = 1 9,jend = 7, clam(l) =. 002, nedge=I, mgsteps=5 S 

The first pass, with equal maximum and minimum spacing, spreads out the / grid lines evenly; the original 
grid points were more densely distributed in the curved section, leaving fewer grid lines in the inlet area. 
This is a good example of using the program to improve an initial grid, with no flow-field adaption 
involved. The second pass steps in the i direction and adapts easily throughout the entire grid. The third, 
fourth, and fifth passes perform the adaption in the j direction. The very different flow features in the 
blunt-body region (blunt-body shock) compared to the features in the inlet region (Mach stem and reflected 
shocks) indicate dividing the adaption domain into these two zones; i < 70 and i>70. Only the blunt-body 
section (i>70) is adapted in pass three, with all default-control parameters. The adaption in the inlet 
domain starts at jst=32 in order to pick up the Mach stem along the lower wall. After stepping through the 
triple-point region (the intersection of the cowl shock and the reflected normal shock), the redistributed 
points do not spread out sufficiently, so the adaption is stopped at line 19 and a fifth pass is started at line 
19 with a decreased value of CLAM. Since this pass starts internally, it is necessary to set MGSTEPS in 
order to merge smoothly from the already adapted region. Note: save=.f has been used on the first four 
passes to prevent excessive file generation. 


Case 5. Subsonic Impinging Jet 

This example is to show the use of the SUB parameter. The original grid, with 231x100 grid 
points, and the corresponding Mach contours are shown in Figs. 21(a) and (b). It was demonstrated that 
adapting the grid enabled the user to decrease the density of grid points, and thus speed up the flow solver. 
Adapting the grid and reducing the number of points can be accomplished in the same execution run. The 
first step is to ensure that the dimension parameters are increased to ID=231, JD=231, KD^l, and 
IMX=23I. (Note that SAGE will inform the user of the correct dimensions, if necessary.) To remove every 
other point from both coordinate directions, the input-control file requires two passes for the adaption: 

Snamel (wod=.t.,sub=l,rdsmax=4.0,rdsmin=. l,clam(l)=.005Jst=I00Jend=l,nedge : =I, indq = 7 $ 
Snamel sub= 1 ,istep=.t.,noup=.t. $ 

Both passes use the SUB parameter, reducing the number of points to a 1 16x51 grid as shown in 
Fig. 20(c). The example shows that the SUB (as with the ADD parameter) can be used in conjunction with 
an adaption (as in the first pass) or simply as a method to reduce the number of points in the original grid 
(as in the second pass). It should be noted that if the two passes had been in reverse order, the input value 
of jst would need to be modified by the user to jst=5L To reduce points in a subset of the adaption domain, 
LSTSUB and LENDSUB should be used; however, in this example the default values of the entire adaption 
domain are appropriate. 
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Figure 21. Subsonic impingement jet. (a) Initial grid; 

(b) initial mach contours; (c) adapted grid, with reduced points. 


Case 6. Axisymmetric Plume Flow 

This concluding 2-D case is presented to indicate the powerful effect of the adaptive grid process on 
the final flow-field solution. Figure 22(a) shows the initial grid and the reflected computed solution (in 
Mach contours) of an axisymmetric nozzle-plume flow described by Venkatapathy and Feiereisen (1989). 
This initial solution has not developed sufficiently to capture the final flow features: the outer shear layer, 
barrel shock, Mach stem, reflected shock, and the triple-point shear layer. Three iterations (through both 
adaption and flow solver) were made to produce the final adapted grid and solution shown in Fig. 22(b). The 
definition of the flow features is greatly improved. Figure 22(c) shows the accuracy of the final solution: 
the lower picture is a shadowgraph of the actual experiment and is almost mirrored by the computed 
solution. 
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Figure 22. Axisymmetric plume flow, (a) Initial grid and flow solution (as Mach contours) showing 
underdeveloped flow features; (b) final adapted grid and flow solution after several iterations; 

(c) comparison of computed solution with experimental shadowgraph . 
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3.2 Three-Dimensional Examples 


The adaption of 3-D grids is more complicated, not only because of the additional torsion control, 
but because of the increased choice of stepping and marching directions. For 2-D adaption, four choices of 
adaption are available: stepping in the / direction, the j direction, or both (in either order). Each of these 
options will produce a different adapted grid. For 3-D adaption, it is also necessary to choose a plane- 
stepping direction and, in theory, to consider up to three possible directional passes. In practice, it has been 
shown that a one-directional pass will frequently provide a sufficiently adapted grid, and that a two- 
directional pass is the maximum required. The following case shows the effect of single- and multi- 
directional passes. 

Case 7. Tandem Fuel Injectors in a Supersonic Combustor (Rockwell Model) 

Figure 23 represents a two-slot, tandem fuel injector arranged behind a backward facing step in a 

supersonic stream that models the 
Rockwell supersonic combustor. The 
fuel injection through the slot nozzle 
creates a complicated three- 
dimensional flow pattern. Since the 
outflow is supersonic, fuel injection 
normal to the main flow produces 
strong shock waves and streamwise 
separation in the vicinity of the slots. 
In addition, the backward facing step 
locally creates a subsonic flow ahead of 
the slots. In this example, grid 
adaption is limited to the region 
downstream of the backward facing 
step. The grid and initial solutions for 
the fuel injection problem were 
provided by J. Wang. 

Example 1: One adaption pass showing 3-D effect. This example is used to demonstrate how a one- 
directional adaption changes the grid-point distribution, not only on the chosen adaption plane, but also on 
the two cross planes. The initial grid (80x31x61) given in Fig. 24(a) shows three selected planes from 
Wang's initial grid, one in each of the coordinate directions, where i-40,j^l and k=30. For clarity, these 
planes are separated (but with their original orientations) and shown in Fig. 24(b). Figure 24(c) presents the 
corresponding Mach contours obtained from the solution on the initial grid by the flow-field code. The 
adaption is performed by redistributing the points on the first j plane and then marching in j planes 
(obtained by setting ikplane=.t .). Within each j plane, stepping is in the i direction (since istep-A.) and the 
adaption is therefore performed along the k lines. The input parameter file used to create the adapted grids 
shown in Fig. 24(d) is: 

$ncimel ikplane = t. , istep = /. , rdsmax=4. 0 , rdsmin =. 2, indq = 7, clam( 1) = 0002 , clam (2) =.0005, nedge = 1 $ 

As requested, the j plane has clearly adapted to the Mach number gradients in the k direction, but so 
has the / plane, since this was the stepping direction within planes. However, the curvature of the planes 
themselves has not changed. The k plane shows a different effect: points have not moved within the plane, 
but only up and down (again, the k direction), thus changing the curvature of the plane. This example 
clearly shows that although the adaption is in only one direction, all planes are affected. It is thus quite 
possible for one pass to adequately adapt the grid for re-input to the flow-field code. 

Example 2: Two-directional pass. This combustor problem is one that benefits from a two-directional 
pass. By looking at the constant j plane in Fig. 24(d), it can be seen that a second adaption to redistribute 
the i direction would be appropriate: this can be effected by adapting the j planes, stepping in k lines. The 
second set of adaption parameters are 

Snamel ikplane =.t., kstep = /. , rdsmax=4 . 0 , rdsmin = 2 , indq = 7, clam (1) =. 0002, clam (2) =. 0005 , nedge = 1 S 



Figure 25. Staged transverse injectors in a supersonic 
combustor, showing computational domain. 
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Figure 24 . Three planes from initial grid l (a) Actual location, i=40,j~ I ,k -30; (h) same planes ; separated 
for clarity ; (c) initial Mach contours; (d) single direction adaption . 

Figure 25(a) shows the result of performing this adaption "on top" of the adaption already shown in 
Fig. 24(d); points have now been redistributed in the / direction as well as the k direction. However, it should 
be noted that the / direction can also be adapted by marching in j within constant k planes by using: 

Snamel if plane t. jstep t. , rdsmax=4. 0, rdsmin =. 2 , indq — 7, clam( 1) =. 0002 y clam (2) =. 0005, nedge= 1 S 

This adapted grid is shown in Fig. 25(b) and shows a very similar redistribution on the j plane, even though 
this is not the adaption plane. 
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Figure 25. Two-directional adaption, using grid shown in Fig. 24(d) as input, (a) Adapting ik planes, 
stepping in k lines; (b) adapting ij planes, stepping in j. 
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Example 3: Torsion parameter between planes, X\ The 2-D examples show how the torsion 
parameter (X) affects the adaption within a plane: it controls the magnitude of the torsion term and the 
larger the value of X, the smoother the resulting grid, but with less grid redistribution. This effect is 
demonstrated in Case I in the 2-D section. Unfortunately, the "base" value of X varies for each problem: 
the default value in the code is .01 but the user may have to change this value by as much as one, two, or 
even more orders of magnitude. X * is the analogous torsion parameter between planes, and acts in a similar 
manner to X, but restricts the movement of points from plane to plane to maintain smoothness in the 
cross-planes. If X* = 0, there is no torsion control between planes, and planes will be adapted as independent 
entities. (In fact, setting X*=0 is one way of handling periodic planes.) In the code, X* is defaulted to 
.0001 and this example illustrates the effect of this parameter. 



Figure 26. Effect of varying X*. (a) Initial grid showing j=l plane; (b) adaption with clam(2)—.001; 

(c) adaption with clam(2)=.0001. 

The same combustor-flow problem from Ex. 1 is used with different input control parameters. This 
time, we are adapting k planes and comparing the effect of X* on the crossing j planes. Figure 26(a) shows 
the initial distribution of two k planes (the upper and lower surfaces) joined by the first j plane. For this 
particular application, the first k plane would normally not be adapted, since the maintenance of its initial 
spacing is a high priority. However, to illustrate the adaption procedure, in this example the first plane is 
adapted. The input parameter file giving the result shown in Fig. 26(b) increases the value of X* to .001 : 
Snamel ijplane-. t.jstep-. t.,clam(l) =. 001, clam(2) =. 001, rdsmax=4. 0, rdsmin =. 2, nedge =1 S 

This input requests marching in k planes (since ijplane=.t.) and stepping in the j direction within the k 
plane. The only parameter with any control of the j plane is X* ( clam(2 )) and the j plane shows some 
redistribution of points, but insufficient to reflect the flow features. Compare this plane to that in Fig. 
26(c), where the adaption was performed with the identical input parameters except for a decrease in X* to 
.0001, i.e., _ 

Snamel ijplane - t. ,jstep=. t.,dam(l) =. 00 1, cl am (2) =.0001, rdsmax=4. 0, rdsmin =. 2, nedge = / $ 

The effect on the smoothing between planes is very evident and indicates the influence of the torsion 
parameter between planes. 


Case 8. NASP 3-D Nozzle Simulation 

Example 1. Two-directional adaption. Figure 27 shows the geometry of a wind tunnel experimental 
model for the National AeroSpace Plane (NASP) called the Single Expansion Ramp Nozzle (SERN). The 
model is inserted into a hypersonic test section with cold air injected at supersonic speed through the 
nozzle. Superimposed on the ramp section is an outline of the computational grid which is detailed in 
Fig. 28(a). This 3-D flow-field grid (41x60x90) defines the nozzle and after-body region of the model 
that is used in the computational experiments. Shown in the figure is the downstream outflow plane at 
i-41 , the lower grid at k=3I (part of which is the upper surface of the after-body ramp) and the symmetry 
plane at j=l . Additional patched grids are not shown, but there is a grid boundary that must be matched 
from j=41 to j=60 along the k=31 plane with the grid that is stored in k=l to k=30, and also a matching 
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i=I plane. Adaptions were performed from plane i=2 to plane i=41, stepping in both j and k directions 
within the plane, using the initial solution shown as Mach contours in Fig. 28(b). This complex case utilizes 
many of the available input options, and the two-pass adaption parameters used to produce the adapted grid 
in Fig. 28(c) were: 

Snamel jkplane=.t.,kstep=J.,rdsmax=2.5 ) rdsmin=.25,kst : =32 t kend=83,indq=7,nedge=l , 
clam(l )= . 1, clam (2) =.01 , mgsteps=4, mgpls =4, march =. t.,ct(l) = . 7 5, save =./ S 

Snamel jkplane =. A jste/? =. t.,rdsmax=l. 25, rdsmin 25, kst =35, indq=7, nedge=2 , wg/ =5, 
clam(l)-. 1 , cl am (2) =. 0 1, mgpls = 4 , ct(l) = 75 5 

The first set of parameters requests the 
adaption of each / plane (jkplane=.t .), 
stepping in the £ direction (kstep=J.) 
within the plane. Mach number ( indq=7 j is 
the adaption parameter. Adaption begins 
on line £=3/ with the merging technique 
invoked ( mgsteps=4 ) within the plane. 
This will retain the k=31 line to match the 
existing lower boundary and ensure that 
the adaption parameters are filtered for 
the next four adaption lines. Both torsion 
parameters and the directional parameter, 
ct(l), are larger than their default values 
since a previous test run showed that too 
much adaption (and thus loss of 
smoothness) occurred using the default 
values. Another parameter used in this example is MARCH. This is invoked from line kend=83 : after this 
line, the flow has no gradient features and is mostly numerical 'noise, 5 and MARCH simply presents a more 
attractive grid. Since the adaption algorithm normalizes the weighting function, 'noise 5 can produce 
unnecessary redistribution. Due to adjacent grids, it is not appropriate to adapt the i~l plane, thus a smooth 
transition is required between the first plane and subsequent planes. This is controlled by mgpls=4 : the 
default start plane of ist=l was not adapted and the next three planes were more gradually adapted than 
they would have been if MGPLS was not requested. Example 2 of this case gives a detailed account of this 
plane-merging process. The SAVE parameter was set to .false, to prevent the output datasets from being 
written. This should be done only when more than one adaption pass is made in the same computer run. 

The second adaption set requests adaption of the same i planes, but stepping in the j direction 
within planes. The first adaption point is k=35. This was chosen to retain the first few points in the 
boundary layer. In addition, the parameters nedge=2 and mgl=8 are used to filter this boundary' layer 
spacing into the next eight grid points on the line. The value of nedge was changed from 1 to 2 for this 
pass because there was no need to maintain any spacing at the outer edge, where no gradients are found: 
these points are better used within the body of the grid. However, the adaption still left an unnecessary 
number of points in this outer region, so a final pass was made: 

Snamel jkplane =. /. jstep - . /. , noup =. /. ,sub= 1, Istsub =82, add-1, lstadd—5 1, lendadd=55 $ 

Here, no adaption is performed ( noup-.t. ) and four lines are moved by using the SUB and ADD parameters. 
The final adapted grid is given in Fig. 28(c). This was input to the flow solver, producing the final solution 
shown in Fig. 28(d). 

ft should be noted that this problem is a multigrid problem and that adaption of grids with matching 
planes can be facilitated by using the export and import capabilities described elsewhere. 



Figure 27. SERN experimental model and the 
computational space in the plume region. 
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Figure 28. NASP nozzle plume flow, (a) Initial grid; (b) initial Mach contours; 

(c) adapted grid; (d) final solution using adapted grid. 

Example 2: Merging from non-adapted planes to adapted planes. The nature of three-dimensional 
problems frequently necessitates the maintenance of boundaries: either for multiple grids or to preserve 
solid geometry' walls and even for maintaining boundary layers. For the NASP nozzle case described here, 
the plane at /=/ must match another plane on the adjacent (not shown) grid. (Note that this plane may 
have received data via the export/import process.) There is also a dense grid spacing around the nozzle exit 
for defining the boundary' layer region that should also be retained. The way to handle these two situations 
is to not adapt the first plane, but to adapt the contiguous planes in a merging manner. This merging is 
needed to create a smooth transition in the cross plane, not in the adaption plane itself (which is handled 
by MGSTEPS). 

Figure 29(a) shows the initial /=/ (and thus the identical i=2 plane) for the 3-D grid shown in 
Fig. 28, and Fig. 29(b) contains the corresponding pressure gradients for the i=2 plane. If the input 
adaption parameter ist=2 is used, with no merging process invoked (i.e., the first plane is ignored) the 
resultant adapted i=2 plane is shown in Fig. 29(c). Although this grid is smooth and nicely adapted to the 
flow, the cross plane, j, is not (compare Fig. 29(c) with Fig. 29(a)). In addition, points have been drastically 
pulled away from the nozzle region. By invoking the MGPLS option, this cross-plane unevenness can be 
dampened. The technique is analogous to that used by MGSTEPS: both A* and C* are modified to increase 
the restraint of movement of points between planes. These values will return to their original input values 
within a certain number of planes, as requested by MGPLS. The following is the input parameter file used to 
create the adapted grid (i=2 plane) shown in Fig. 29(d): 

Snamel jkplane =. /. .jstep =. /. , rdsmax=J.25, rdsrnin =. 25, kst=3 5, indq=6, nedge =2, mg l =8, 
clam(l)-. l,clam(2) =.01, ct( 1) =. 75, mgpls=4 $ 
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The difference between Figs. 29(c) and (d) is striking: the boundary layer spacing is now maintained and the 
adaption process is just beginning. Although mgpls=4 in this example, setting mgpls = I will also have 
considerable effect. By setting MGPLS at all, a request is made to not adapt the first requested plane (in this 
case, m 7=7, by default) but to use the grid on that plane as a control for the subsequent planes. Thus, 
mgph=l, is(=l is_not the same as setting ist=2. 





(d) 


Figure 29. Merging planes using mgpls . (a) Initial grid at i—1 (and i=2); (b) initial pressure contours; 
(c) adapted grid at i=2 with ist=2; (d) adapted grid at i=2 with ist=l and mgpls =4. 


Case 9. Aeroassist Flight Experiment (AFE) Vehicle 

This example is the hypersonic, non-equilibrium flow around the forebody of the aeroassisted flight 
experiment (AFE) vehicle. In Fig. 30(a) the initial grid configuration (35x23x49) is shown around the 
body in the form of j planes 2 and 22, and the outflow plane at i=35. The Mach contours seen in Fig. 30(b) 
were computed by the flow-field solver of Palmer (1990) using the non-adapted grid. This grid was then 
adapted with respect to the Mach contours, and the resulting grid is seen in Fig. 30(c). The redistribution of 
points within the blunt-body shock region is clearly shown. Due to the smooth shape of the shock, 
adaption was performed only in one direction: marching in j planes and stepping in the i direction within 
planes. The input data set used to create the adapted grid in Fig. 30(c) was 

Snamel lnsing=l f ikplane t. , istep=. t. , indq = 7Jst=2Jend=22, rdsmax = 4 . 0, rdsmin =. 1, 
kst=15, orthe=f,nedge=l f clam( I) =. 001 $ 

This is the first example to use the LNSING parameter. It is set since all j planes emanate from a 
single line at i=I. (This forebody singular line is used in grid mapping.) However, the dataset still contains 
the grid points for line i~l as if it were a separate line in each plane. lnsing=l ensures that the first line in 
the first plane is adapted, and that all /=/ lines on subsequent planes are set identical to this adapted line. 
The parameter orthe was set to .false. , which removes the orthogonality constraint at the outflow line of 
i=35. The default value of .true, would have created a false turning of the adaption at this location. This 
adapted grid was then re-input to Palmer's flow-solver which showed a considerably sharper blunt-body 
shock feature (1990). 
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Figure 30. Hypersonic forebody flow (AFE). (a) Initial grid around the forebody; 

(b) initial Mach contours; (c) adapted grid. 

3.3 Multiple-Grid Examples 

Section 2,4.3 describes some of the complexities of multiple grids. The first case here reproduces 
the stylized example shown in Fig. 1 1 to illustrate the transfer of planes from one grid to another. The 
second case is a generic NASP configuration (Davies and Deiwert, 1993). 

Case 10. A Simple 3-D Grid 

Example 1: Single plane transfer. Figure 31 is a reproduction of Fig. 11 from Section 1. This example 
shows a matching common plane that is stored separately in each grid, but contains identical data. We will 
define the size of each grid as (ij,k) of (40x15x30). If the grids had been separated as single grids, as 
shown in Fig. 31(b), the only choice would be to start the adaption process at the common plane and to 
march in opposite directions. This would provide identical adaption points for the common plane. For 
some problems this could be quite acceptable, but marching in the same direction more often produces the 
preferred result. 



To adapt multizoned grids as shown in Fig. 31(c), a feature is available to transfer data from one grid to 
another. The adaption process in Fig. 31(c) is to: 
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(1) adapt grid 1 up to the common plane (j=15)\ 

(2) transfer data from the last plane in grid 1 to the first plane in grid 2; and 

(3) adapt grid 2, using the merging planes option (i.e., mgpls nonzero) for the first plane. 

This sequence is described by the following input control list: 

Snamel mgr id = 1, ikplane =. /. S 
Snamel export=.t.,mgrid=l , ikplane-. t.,mplane=l 5 S 
Snamel import=.t.,mgrid=2,ikplane=.t.,mplane=l S 
Snamel mgrid=2, mgpls =2, ikplane =. t. $ 

If only a subset of plane 15, grid 1 was to be transferred to grid 2 (as shown in Fig. 32), then the 
range would be described by entering values for parameters is,ie,ks,ke (i.e., those relating to the ikplane) for 
both grids by 

Snamel export =. /. , mgrid= 1 , ikplane-. L, mplane =15, is =30 S 
Snamel import = /. , mgrid=2 ikplane =. /. , mplane =l,ie = 10 $ 

In this case, points 30</<40 and all k points from the 15th ikplane in grid 1 are transferred to points 
1 < / < 10 in the first plane of grid 2. (Note that only is=30 on the export list and ie=10 on the import list 
are required since the other range values are the same as the default values.) 

Note for 2-D files . The same concept holds for 2-D files but care 
must be taken since only a line is being transferred. Because only one 
plane exists, ijplane=.t. and mplane=l . Only is, iejsje can be used and 
one of these two sets must have equal first and last values. Finally, if 
the export card is the first card in the adaption set, include twod=J. 

Example 2: Multiple transfers* It is possible to handle more than 
one transfer within the same run of the code. In the example shown in 
Fig. 32, it is probable that the lighter shaded domain of plane 15, 
grid 1, matches to a third grid, not shown. For this example, the two 
sets of transfers can be entered sequentially: 

Snamel export = t. , mgrid=l, ikplane =. /. , mplane = 15, is -30 S 
Snamel import =. /. , m grid =2 , ikplane — . (. , mplane - l,ie=10 S 
Snamel export =. t. , mgrid= 1, ikplane =.t., mplane = 1 5, ie =30 $ 

Snamel import =. t . , mgrid= 3 , . ikplane =. t. , mplane =1 S 

Case 11. Generic NASP Configuration 

This case is the result of NASP Government Work Package #20, by Davies and Deiwert (1993). A 
two-part multiple grid defines a generic NASP vehicle configuration: one grid (35x59x51) for the nose 
region and one grid (also 35x59x51) encompassing the remainder of the body. Figure 33(a) shows the 
vehicle surface and three cross planes (Jk planes), and Fig. 33(b) shows the corresponding Mach contours 
that were used as the adaption function. The i direction is through the vehicle centerline, j steps away from 
the body and k marches around the body from the lower surface. Adaption begins at the nose and marches 
plane by plane towards the end of the body, transferring data at the matching plane (grid l , jkplane 35 to 
grid 2, jkplane 1), The input control file used to produce the adapted grid shown in Fig. 33(c) is 

Snamel mgrid= 1 Jkplane=.t,,ist=2,kstep=.t.,kst=5 1 ,kend= 1 Jst= 1 5,indq=7,nedge=l , 
rdsmin =. 1, rdsmax=3. 0, clam(l ) =. 001, ct(2) 8, orthe(l )=.f, orthe(2) =./ $ 

Snamel export=.t.,mgrid=lJkplane=.t.,mplane=35 S 
Snamel import=J.,mgrid=2jkplane=.t.,mplane=l $ 

Snamel mgr id=2, jkplane =. t. , kstep=. t. , kst =5 1, kend= 1, indq = 7, mgpls =4, nedge = 1 , 

rdsmin =. 1, rdsmax-3. 0, clam (1)=.00I, ct(2) =. 8, orths (2) =.f, or the (1) =.f, or the (2) =.f S 

Some of these adaption parameters need explaining: ist—2 is the starting plane since the nose at ist=l is a 
plane collapsed to a single point; kst=51 and kend=l produces the best adaption direction for kstep 
(compared to the result using the defaults kst=l and kend=51 ); jst=15 preserves the boundary layer; and 
ct(2)=.8 (more straightness for the between-pianes torsion parameter) produces a smoother adaption when 
looking at the other planar directions. The export and import processes are straightforward. For the 



Figure 32. Matching subsets . 
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adaption of grid 2, mgpls=4 and orths(2) = / both contribute to maintaining continuity across the matching 
plane. 



(c) 


Figure 33. Generic NASP vehicle configuration, (a) Initial grid; (b) initial Mach contours; 

(c) adapted grid. 

3.4 Recluster and Boundary Movement Examples 

Reclustering is an alternative to adaption. It was initially added to give the user a second option for 
point redistribution after the outer boundary was moved. (The first option was to proportion the new grid 
spacing to the old grid spacing.) However, the same reclustering algorithm proved to be a useful tool to 
redistribute the grid points without any boundaiy movement. An example would be to convert a grid 
developed for an inviscid calculation into one appropriate for viscous flows. The first example below 
demonstrates the reclustering technique without boundary movement; further examples demonstrate 
boundary movement options. 


Case 12. Space Shuttle 

A single grid around the Space Shuttle, shown in Fig. 34(a), is used for these examples. The i 
direction runs from the nose to the tail, and the k direction is from the body surface to the outer boundary, 
thus the mesh seen in Fig. 34 is an ikplane. This initial grid is equally spaced to more clearly demonstrate 
the reclustering process. Note that SAGE can transform any grid into an equally spaced grid by using 
rdstnax=1.0, rdsmin=1.0, clam(l)-.0, clam(2)-.0. 

Example 1. Reclustering entire domain with no boundary movement. To create a grid with the 
stretched spacing shown in Fig. 34b, the following input dataset was used: 

$namel lnsing= I . ikplane =. /. , islep=. t. , r edits t= 1, mvbound=0, noq =. t . , dsw 0001 ,S 
The parameter Insing is set because of the singular line at the nose, and mvbound=0 signifies no boundary 
movement. No Q file is needed to recluster the grids, so noq is set to true (however, this implies that any 
existing Q file will not be interpolated onto the new grid points). Setting reclust=l is the same as 
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reclust=kmax (ikplane and istep determine which lines we are reclustering, in this case, lines of constant k\ 
hence all the points will be reclustered. The value of dsw (.0001) specifies the actual size of the first spacing off 
the wall, not a normalized value. 



Figure 34. Space Shuttle grid, (a) Initial, equally spaced grid; (b) new grid, reclust~l; 
(c) new grid with reclust=l and dse=1.0. 


Figures 34(b) and 34(c) show the effect of the reclustering. Points are tightly spaced next to the vehicle 
surface for both cases. In Fig. 34(b), the outer Vinokur parameter, dse, contains the default value of 5.0 (which 
sets the outer spacing to five times the average spacing); however Fig. 34(c) was created by inputting dse= 1.0, 
which retains a more dense spacing at the outer boundary. 

This grid is an example where the outer boundary is close to the vehicle surface at the nose and then 
rapidly expands away from the vehicle surface downstream. In this situation, the absolute size of the first wall 
spacing, As i=] , at the nose may not be appropriate for A s jzzimilx . In addition, such tight spacing at the nose can 
create difficult conditions for the Vinokur solver further along the body, producing an unsmooth grid. This is 
clearly shown in Fig. 34(d) where the value of .00001 has been used for the wall spacing all along the body. 



Figure 34 concluded, (d) Constant wall spacing, dsw-. 00001; (e) expanding wall spacing . dsw =999, 

dswl~ 00001, dsw2=. 0001. 

To create the smoother grid seen in Fig. 34(e), the user can set dsw=999 and input two values for the 
wall spacing, one at the nose (dsw I =.00001) and a larger value at the tail (dsw2 =.0001), by using the input 
control file: 

Snamel lnsing=l j ikplane=.t.,istep=.t. i reclust=l,mvbound=0,noq=.t., 
dsw =.999., dsw / = . 0000 1, dsw 2 = 000 1 $ 

If reclustering is required within a subset of the grid, as for example, within the boundary layer, the 
input control file would be the same, except for reclust , which would now be set equal to the number of points 
to recluster, where reclust<kmax. 
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Example 2: Percentage boundary movement. To move the outer boundary by a percentage of its original 
location, mvbound is set equal to the percentage required. Figure 35(a) shows the same original grid as Fig. 
34(a). The input file 

Snamel lnsing=l,ikplane=.t.,istep=.t.,reclust=0,mvbound=-20,noq=.t., S 


produces the grid shown in Fig. 35(b): the negative sign on mvbound pulls in the outer boundary by 20%. 





Figure 35. Percentage move in outer boundary, (a) Original grid; (b) mvbound=-20; 

(c) mvbound=—20, nsm-0. 

If mvbound=20 had been used, then the boundary would have moved out by 20%. This example also 
provides a glimpse into the smoothing parameter nsm, which has a default value of 10: the smoothing has 
caused a blip (shown by the arrow) at the change of boundary length over the shuttle window. The case was 
rerun, adding nsm=0 to the same input file above, giving the result seen in Fig. 35(c). In both cases, reclust=0 
provides a new grid distribution proportional to the original spacing. 

Example 3: Outer boundary movement parallel to shock. This example demonstrates moving the outer 
boundary to match a shock. Figure 36(a) is the initial grid and the computed density contours are shown in Fig. 
36(b). The following input dataset was used to produce the grid shown in Fig. 36(c): 

Snamel lnsing= l ,ikplane=.t.,istep=.t.,reclust=l,dsw=.001,mvbound=999,indq— 1 $ 

When mvbound=999, the code is told to use the q variable of record (in this case indq=l, density) to establish 
the new location of the outer boundary. The code will search from the outer edge until a jump in density 
gradient is found. In the example, the boundary is moved inwards, but it is possible for a boundary to move 
outwards if the shock crosses through the outer edge. However, if the outward movement is not satisfactory, one 
technique is to move the boundary out by mvbound=+% and then re-run with a mvbound =999. This will ensure 
that the entire shock is captured downstream and points are not wasted in the nose region. 



Figure 36. Outer boundary > movement, (a) Initial grid; (b) density contours; 
(c) new grid, with outer boundary moved inwards and points reclustered. 
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There is another option to move the boundary to match a contour line. If mvhound > 1000.0 , then the 
boundary will run parallel to a particular contour line whose value is mvbound- 1000.0. Choosing the variable 
(i.e., indq) and the value of the variable needs care. For example, consider the following input: 

$namd Insing^lJkplane^.t.Jstep^.t. ,reclust^l,dsw=. 001 ,mvbound= 1400. 0,indq=5 5 

Note that the fifth Q variable in this example contains temperature. The temperature contours and the 
resulting grid are shown in Figs. 36(d) and (e). The value of mvbound= 1400.0 requests the outer boundary to 
run parallel and close to the 400° contour line. Care must be taken to ensure that the chosen value exists along 
all lines to prevent an unacceptable result. 



Figure 36 concluded, (d) Temperature contours; (e) outer boundary line follows the T = 400 ° contour. 
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5. SAGEv2B, THE BLANKING VERSION 


This section describes the blanking version of SAGE. SAGEv2 was modified to take into account 
the BLANK feature available in PLOT3D, and was renamed SAGEv2B. The current SAGEv3 described in 
this document does not accept the input parameter BLANK described below. Users interested in this option 
should contact the authors for SAGEv2B. When time permits and/or if there is sufficient demand, this 
feature will be added to the latest version. However, the following description will remain valid. 

5.1 Algorithm for Blanked Overset Grids 

As mentioned in Section 1.10.1, an overset blanking grid system can produce complicated grid 
structures. Grids in the same system can have different functions. For example, a fine, detailed grid may be 
created around each physical object, along with a larger, less dense grid that overlays the entire domain. In 
fact, several layers of nested grids may be appropriate. Nevertheless, as with any grid generation method, 
grids are still created around the physical objects with no a priori knowledge of the solution. Since the 
accuracy of the flow solution is dependent on the grid, the need for an adaptive grid scheme is still 
important. 

One obvious option is to ignore the blanking flag and adapt the grid using whatever data is available 
in the blanked regions. However, in most cases these values are not relevant to the solution. Therefore, a 
change in the algorithm and some careful coding is required. 

The basic adaption equation is given by Eq. (7): 

co, As, - CO,., As,., + /( x, ) + f(V,) = 0 

where 5 is the streamwise length along the adaption line, oo is the tension force between points on a line, x 
is the within-plane torsion force, and \|/ is the between-plane torsion force. The X torsion force is a 
function of the torsion vector defined in Eq. (16), where the constants /„ and C, may change from line to 
line but are constant along a line, and hence are defined as single variables. When blanked regions are 
involved, there will be active points along a j line that may be blanked (i.e., inactive) on the j-1 line. For 
the straightness vector e , data is also required from the j— 2 line. A similar situation exists for the between- 
plane torsion vector, t* . A simple method was chosen to handle the computation of the torsion vectors 
when blanked points are found. First, the values of C, and t n were made dependent on /. If the point (i.j.k) 
exists, but (ij-l,k) is blanked, then t n =1 and C, = 0 will place all the torsion force onto the orthogonal 

vector b, acting on the j line. If the points ( i,j,k j and (ij-l,k) exist, but (i,j-2,k) is blanked, then the 
straightness vector is amended to include the j—1 point only. The same changes are made to the 
coefficients of the vector, depending on the existence of the (i,j,k-l) and ( i,j,k-2 ) points. 

Another significant difference is the number of points on a line. With the blanking feature invoked, 
the number of / points along a j line can differ line by line, varying from the user’s input domain. The user 
specifies the maximum value of (istjend), (jstjend), and (kst.kend) for the adaption domain; the code will 
recompute the start and end points of the adaption line for each step, within the user’s specified limits. 
When blanking occurs in the middle of a line, the two segments must be adapted separately. 

Special handling occurs at boundaries. For edge boundaries, continuity is maintained between 
adaption and physical boundaries. A question occurs along a boundary created by a blanked region: is this 
edge to be externally influenced or not? It is currently assumed that edge control will not occur at a blanked 
boundary unless the blanked boundary point is closer to the outer domain edge than the value of the 
corresponding MG value. For marching boundaries in blanked applications, the code determines the start 
and end points along the last adapted line (either jend.k or j,kend, depending on MARCH or MARCHPL). If 
a line to be interpolated starts (and/or ends) within the range of the last line, then only the appropriate 
points will be used from the last line. If there are more points along the interpolation line, only the 
available range from the last line will be used: the remaining points will retain their initial locations. 

The addition or subtraction of points has been handled very simply. On request, points will be added 
or removed over the entire adaption domain, regardless of blanked regions. The value placed in the 
blanking array will reflect the proper position of each point, indicating whether it is in the blanked region 
or not. 
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The proportionality coefficient, A, controls the magnitude of the torsion parameter. When a 
situation occurs that drastically increases the number of points in a line, the value of X, although 
appropriate for the overall adaption, may permit too much movement when integrating the blanked region 
into a non-blanked region. For this reason, X is locally magnified to prevent excessive point movement. 

5.2 Execution of SAGEv2B 

SAGEv2B is the only version of SAGE that recognizes the blanking option that is found in 
PLOT3D grid files as an additional variable, called IBL. The changes to the adaption algorithm to handle 
blanked areas are discussed above. (Note: Tf IBL is defined in the PLOT3D input grid file, SAGEv3 can be 
used but the blanking will be ignored and the output grid file will not contain IBL.) 

Invoking the blanking option is extremely easy: the variable bktnk=.t. should be added to the input 
parameter file. However, defining the adaption domain and the marching directions requires much more 
care. With the blanking feature invoked, the number of z points along a j line can differ line by line, varying 
from the user’s input domain. The user must specify the maximum value of ( ist.iend '), (jstjend), and 
(kst,kend) for the adaption domain; the code will re-compute the start and end points of the adaption line 
for each step, within the user’s specified limits. If the blanked region is in the middle of the grid, then 
another strategy must be developed. One option is to divide the adaption domain into two. Unfortunately 
this fixes one of the lines or creates a discontinuous grid. Performing two adaptions, sweeping in opposite 
directions, will provide a better result. The examples section will clarify this. 

5.2.1 User Input Parameter 

There is only one new input parameter: 

BLANK (FALSE) set to .true, to read PLOT3D IBL variable 

BLANK The grid file in PLOT3D has an additional optional variable. Following the X,Y,Z arrays 
is an array called IBL. IBL is the same dimension as the grid and contains an integer value that describes the 
status of the grid point. The standard nomenclature of IBL=1 implies the grid point is active, and IBL=0 
signifies an inactive point. Some users may put other values in TBL, for example IBL=2 to represent a 
border point. This code assumes that if IBL=0, the point is not considered; otherwise, it is assumed to be 
part of the grid and will be counted as an active point on a line. The input values of 1ST, IEND, etc., that 
describe the active adaption domain must contain the maximum required domain; the code will determine 
the active domain for each line based on the values of IBL and these maxima. 

5.3 Subroutines and Variables 

Several previously defined subroutines are amended. For example, all the I/O routines are changed to 
handle the IBL variable (including reading, writing, swapping etc.) Also, variables such as NIPTS, MG I, etc., 
that were constants but now vary are renamed so they can be retrieved, and variables that become arrays 
(such at CT) are defined. Because of its expanded size, the routine INITIAL is split into two: INITIAL and 
GETDATA. In addition, three new subroutines w'ere created. 

BL TST(IB,JB, KB, BLANKIT) 

This routine tests to see if any points are blanked around an existing point ( ib,jb,kb ), for example, a 
‘hanging’ point. If so, this point is also treated as blanked since there is no connectivity for adaption. The 
variable BLANKIT is therefore set to true. 

LNPTS(J,K) 

This routine computes the number of active points on a line for adaption. Blanked points are not used. 
Also the edge merge parameters (MG1 and MG2) are amended if blanking occurs within the line edge. 

TORBL(J.K) 

This routine amends the torsion coefficients C, and t n for (i,j,k) to reflect the status of adjacent points. If 
the point at ( i,j-l,k ) is blanked, then within-line coefficients are changed to C, =0 and r„ =1. If (ij,k-l) 

is blanked, then the within-plane coefficients are changed to C* =0 and f* = I. 
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New variables are: 


AIJ 

BLANK 

BLANKIT 

IB,JB,KB 

IBL 

ISTFJENDF 
MG 1 F,MG2F 
NIPTSJ1 


/COM25/TORCOF/ used to amend X near blanked points 
/COM25/input/ invokes blanking option 

/argument/BLTST/ sends flag to LNPTS indicating status of current point 
/arguments/BLTST/ location of current point from LNPTS to BLTST 
/COM25/input/ array containing blank flags on input grid file 
/COM25/GETDATA/ saved input values of 1ST and IEND 
/COM25/GETDATA/ saved input values of MG1 and MG2 
/COM25/LNPTS/ saved value of NIPTS on j-1 line 


5.4 Example Case: Access-to-Space 


The blanking version of the SAGE code was developed to assist the Access-to-Space project, 
particularly the plume and base-flow computations. This example case consists of a multiple grid with a 
total of five overlaid grids enveloping the base and plume flow regions. Grid number three overlays the 
entire domain and contains many of the features that the new blanking version of the SAGE code must 
accommodate. It was therefore chosen as the test case of the adaption process. Figure 37 shows the 
constant j—2 plane and two constant i planes, 1=48 and i=20, from grid #3 that is dimensioned 
(48x31x61). The arrows in the figure indicate the coordinate directions. 



Figure 37. Grid 3 domain. 


Adaption of a single i plane . 

To examine the methodology of adapting around blanked regions, we first look at a single plane. 
Figures 38(a) and 38(b) show the initial grid and the interim Mach contours calculated on the constant j=2 
plane seen in Fig. 37. 



Figure 38. Plane j=2. (a) Initial grid; (b) initial Mach contours. 
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Since this was an interim solution, many of the flow features are not clearly developed. In addition, 
only 'noise 5 can be seen at the onset of the flow (/ planes 1 through 16). As a result, the first adaption 
procedure was initiated at i=16. This example is complex because it contains a blanked region at a 
boundary, an inner blanked region, merging from two blanking areas to a full grid, and even some 'hanging 5 
points/lines (points that are only singly attached to the grid) that can be seen halfway along the blanked 
boundaries. 

The plane was adapted using the following user-parameter input file: 

Snamel mgrid=3, blanks. t. , ik plane =. t.,istep=. t., jst =2Jend=2, ist= 1 6, i?idq=7 J 
rdsmax-2.5,rdsmin=.4, clam(l)=.05, nedge = 1, mgsteps = 4 $ 

The mgrid=3 parameter informs the SAGE code that the grid and solution files are in multiple grid format 
and that these adaption parameters apply to the third grid. The parameter blanks. true, ensures that the 
code reads the blank (IBL) array and tests for any blanked regions during the adaption of each line. The 
ikplane and istep parameters define the line-stepping and plane-stepping directions. The adaption domain is 
the j=2 plane, / points range from 16 to 48, and all the k points are included (if the domain range is not 
input, physical boundaries are used as default). Indq =7 indicates that Mach number is the adaption variable; 
rdsmax and rdsmin define the relative ma xim um and minimum me sh spacing: clam(l)=.05 increases the 
spring force, cd, to produce a smoother adaption; nedge=l invokes the edge spacing control and mgsteps=4 
invokes the merging algorithm at i=16 to moderate the discontinuity between the adapted and non-adapted 
region. All of these parameters are described in detail earlier in this user guide. 



Figure 39. Adaption from left to right , starting at i=16. 

This was a single-sweep adaption and the resulting grid is shown in Fig. 39. It can be seen that the 
upper left side of the grid has not been adapted. This is because the adaption algorithm can only handle one 
segment of a line, and since the line point count began at k=l , only the segment between the two holes was 
captured. In order to capture the upper segment, a second sweep is required, this time defining the adaption 
domain with kst=6 7, thus ensuring the line point-count algorithm starts looking for active points at k=61 . 
The following two-sweep input parameter file produced the result shown in Fig. 40. 

Snamel mgrid=3, blank .t., ikplane -.t.,istep=.t. f indq=7,jst=2Jend=2, ist= 1 6,kst=6 1 , 
kend=I,rdsmax=2.5,rdsmin=.4, clam(l)=.05, nedge=l, mgsteps-4 S 
Snamel mgrid=3 , blank— A., ikplane =. istep=. t. , indq-7Jst=2Jend=2 , ist=16, 
rdsmax=2.5,rdsmin—.4 , clam(l)=.05> nedge=l t mgsteps=4 S 

This adaption provides a quite satisfactory result: all points have been adapted to align with the flow 
features and no distortion of the blanked areas has occurred. 



Figure 40. Two-sweep adaption . 
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Adaption of all planes. 

Finally, Fig. 41 shows the complete adaption for all j planes. The input data set is almost 
unchanged: the single-plane restriction on j has been removed and a value for the torsion between planes 
has been input, giving 

Snamel mgrid=3, blank=.t., ikplane~.t.,istep=.t., indq=7,ist= 16, kst=61,kend=l, 

rdsmax =2.5, rdsm in =. 4, clam(l)=.Q5, clam(2)=.05 nedge=l, mgsteps=4, $ 

$name mgrid=3, blank=.t., ikplane=.t.,istep=.t., indq = 7, ist =16, rdsmax =2. 5, rdsm in =. 4 
clam(l)=.05, clam(2)=.05, nedge=I, mgsteps=4, S 

The figure shows the i—20 and i=48 planes and the j=2 and jend planes. All the flow features have been 
captured by the adaption. 



Figure 41. Full adaption of computational grid. 
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